{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Logistic回归与Softmax回归"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 数据准备"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "import itertools\n",
    "import functools\n",
    "\n",
    "class PolynomialFeature(object):\n",
    "    \"\"\"Polynomial features\n",
    "\n",
    "    transforms input array with polynomial features\n",
    "\n",
    "    Example 1\n",
    "    =========\n",
    "    x = [a, b]\n",
    "\n",
    "    y = PolynomialFeature(degree=5).transform(x)\n",
    "    y = [[1, a, a^2, a^3, a^4, a^5],\n",
    "         [1, b, b^2, b^3, b^4, b^5]]\n",
    "\n",
    "    Example 1\n",
    "    =========\n",
    "    x = [[a, b],\n",
    "         [c, d]]\n",
    "\n",
    "    y = PolynomialFeature(degree=2).transform(x)\n",
    "    y = [[1, a, b, a^2, a*b, b^2],\n",
    "         [1, c, d, c^2, c*d, d^2]]\n",
    "    \"\"\"\n",
    "    def __init__(self, degree=2):\n",
    "        \"\"\"Construct polynomial features\n",
    "        \n",
    "        Params\n",
    "        ======\n",
    "        degree: int\n",
    "          degree of polynomial\n",
    "        \"\"\"\n",
    "        if not isinstance(degree, int):\n",
    "            raise ValueError(\"degree should be int type\")\n",
    "        self.degree = degree\n",
    "\n",
    "    def transform(self, x):\n",
    "        \"\"\"Transforms input array with polynomial features\n",
    "\n",
    "        Params\n",
    "        ======\n",
    "        x: ndarray with shape (N, D)\n",
    "           input array\n",
    "\n",
    "        Returns\n",
    "        =======\n",
    "        output: ndarray with shape (N, K+1)\n",
    "           polynomial features\n",
    "        \"\"\"\n",
    "        if x.ndim == 1:\n",
    "            x = x[:, None]\n",
    "        x_t = x.transpose()\n",
    "        features = [np.ones(len(x))]\n",
    "        for degree in range(1, self.degree+1):\n",
    "            for items in itertools.combinations_with_replacement(x_t, degree):\n",
    "                features.append(functools.reduce(lambda x, y: x * y, items))\n",
    "        return np.array(features).transpose()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_toy_data(add_outliers=False, add_class=False):\n",
    "    x0 = np.random.normal(size=50).reshape(-1, 2) - 1\n",
    "    x1 = np.random.normal(size=50).reshape(-1, 2) + 1\n",
    "    if add_outliers:\n",
    "        x2 = np.random.normal(size=10).reshape(-1, 2) + np.array([6., 3.])\n",
    "        return np.concatenate([x0, x1, x2]), np.concatenate([np.zeros(25), np.ones(30)]).astype(np.int)\n",
    "    if add_class:\n",
    "        x2 = np.random.normal(size=50).reshape(-1, 2) + 3.\n",
    "        return np.concatenate([x0, x1, x2]), np.concatenate([np.zeros(25), np.ones(25), np.zeros(25)+2]).astype(np.int)\n",
    "    return np.concatenate([x0, x1]), np.concatenate([np.zeros(25), np.ones(25)]).astype(np.int)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x11b71a790>"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAD4CAYAAADxeG0DAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dd5hU5dnH8e89fXaXXdqCSBMRUYN9JWCPLWDBqK+VqIkFO7Yk5hU1mljx1cRuUFFjVIzYC/ZCLCBLEVFQsFEVpG+Zfr9/zIrs7qwsO+XM7Nyf69oL9szseX5DuefMc54iqooxxpjC5XI6gDHGmPRYITfGmAJnhdwYYwqcFXJjjClwVsiNMabAeZxotGvXrrrVVls50bQxxhSs6dOn/6CqlU2PO1LIt9pqK6qrq51o2hhjCpaIfJvquHWtGGNMgbNCbowxBc4KuTHGFDgr5MYYU+CskBtjTIFzZNSKMU5QTaB1E6H+X5CogcDBSNk5iKuz09GMSYsVclM0dN1VUP8CUJ88UPcYGnoVur6MuMoczWZMOqxrxRQFjS+F+ufYUMQBiEJiDVo/0alYxmSEFXJTHKJzQLwpHghB5MOcxzEmk6yQm+Lg6g4kUjzgAXevXKcxJqOskJvi4N0J3D1pflvIg5T81olExmSMFXJTFEQE6fQQeHcDfEAAXN2RTncjnn4OpzMmPTZqxRQNcVciXf6NxleC1oK7NyLidCxj0pZ2IReRADAZ8Decb6Kq/iXd8xqTLeLuAnRxOoYxGZOJK/IwcICq1oiIF3hPRCap6pQMnNsYY8wmpF3IVVWBmoZvvQ1fmu55jTHGtE5GbnaKiFtEZgHLgddVdWqK54wSkWoRqV6xYkUmmjXGGEOGCrmqxlV1F6AXMFhEBqV4zjhVrVLVqsrKZjsVGWOMaaOMDj9U1TXAO8CwTJ7XGGNMy9Iu5CJSKSIdG34fBA4C5qV7XmOMMa2TiVErPYCHRcRN8o3hP6r6YgbOa4wxphUyMWplNrBrBrIYY4xpA5uib4wxBc6m6BuTh1SjEHodDb8Dri5IybGIZ2unY5k8ZYXcmDyjGkFXnQyxz0HrAA9a9yhacROu4HCn45k8ZF0rxuQZrXsaovMaijhADAjBustRDTsZzeQpK+TG5JvQizTeku5HAtGPc53GFAAr5MbkGylp4YEEyUVGjWnMCrkxeUZKTgSCKR4oB++OOc9j8p8VcmPyjX9/KDkJ8CevzqUUpBPS6T5E7L+sac5GrRiTZ0QEKb8MLf0tRD4C6Qj+vUku929Mc1bIjclT4u4JwaOcjmEKgH1OM8aYAmdX5MaYdkM1npxIhQs8A4tmc20r5MZkmSbqIPwm6Drw7Yl4+jkdqV3SyDR0zWjQhjH4UgGd7ka8v3A2WA5YITcmizQyA119BqCgcUDRkuORDmOK5moxFzSxCl195kazYQGtQ1edApX/RVwtjc1vH6yP3JgsUY2iq88GrQGtBUJAGOqehMi7TsdrX+pfaHijbCoO4ddzHifXrJAbky2RGSTXSWmqHq17Mtdp2jVNrABSrEOjUUj8kPM8uWaF3JisSVXEG2g0dzGKgPiGtLC0gQd8g3OeJ9eskBuTLb7dSK6P0oQEkeCInMdp13x7gmcnGi9tEAT/vkgRLGtgNzuNyRKRIFo+Ftb+AYgD0eRVo3cPCNi64ptDNQKxeSCliKd/s8dFXND5/mSXVegZwIMEj4PgkbkP6wAr5MZkkSt4COqbhNY9A7oW8e8PvqG2ZspmSNS/AuvGkBz5E0PdvZBO9yKePo2eJ+JDSkdC6UhngjrICrkxWSbunkiH852OUZA0Oh/W/onkiJ8G8a/QVadC5Zv2htgg7T8FEektIm+LyFwR+VRELsxEMGOM0brHgEiTownQNRCd7kSkvJSJK/IYcKmqzhCRDsB0EXldVT/LwLmNaTWNfobW3JWcou3ZFik7ryhm9bVrie9IecMYgcTKXKfJW2lfkavqMlWd0fD79cBcoGe65zVmc2ikGl15AoTfgPhCCL+JrjwRjXzkdDSTDt/+pNxkQyPg3TXXafJWRjuYRGQrYFdgaorHRolItYhUr1ixIpPNGoOuu5ZkP6r+eAQINRw3hUpKjgT3ljTa4k6CUHIK4u7uWK58k7GbnSJSBjwFXKSq65o+rqrjgHEAVVVV2vRxY9ISm9ficVW1dU0KlEgAujyJ1j0KoUkgHZDS34L/EKej5ZWMFHJJbl3yFPCoqj6diXMas1mkAnR1yuNWxAubuMqQsrOg7Cyno+StTIxaEeABYK6q3pp+JGPaoPT3NO9LDUDp7xwIY0xuZaKPfC/gZOAAEZnV8HVoBs5rTKtJ6SgoOZ7khsWlyV9LjkVKz3Y6mjFZl3bXiqq+B9hnV+MoERdSfjlaNhriS8G9JeIqczqWMTlhMztNuyKuMnBt63QMY3LKCrkxxmSRagStHQ/1E5ObXwSPQErPQlylGWvDCrkxxmSJqia3oIvMZMN6MbXj0fDb0OUZRDJTgm3FGWOMyZboTIh+TKNFv4hAfDGE38pYM1bIjTEmW6Ifg6bYKUpr0cjMjDVjXSsFJJFI4HLZe2++09hitPYeiEwFdw+kdBTi38fpWMYJ7h4g3uTaMI0EEHfmlqSyqlAA3prwHif1PZtfe47nuB5n8MK9r6JqqxzkI40tRlceCfVPJxfvikxFV59Pou4Jp6MZJ/gPSK4N07TUigeCR2SsGSvkee6/T0/l1jPuYcWi5JKdq79fyz//8AjP3/2Kw8lMKlp7F2gdya3dflQP68cmtyszRUXEh3R+HDzbAT7AD+6tkc7/RlwVGWvHulby3INjHiNc17gAhOvC/OuaJxlx7jBbRyTfRKbSuIj/KA7xRahrC7T2Xqh/AXBB8CikbBQi/hQ/Y9oD8fRFuj6LxlcAccS9RcbbsEKe577/NvWSvzWraoiGo/gCvhwnMj/L1T05IqEpjaFSAatOgtiXbNj1pnYcGnkfOj+e8k1ZY18mN8uIzARPb6T0HMQ/tNVxNL4EQm+BuMF/EOLu1sYXZtIl7sqsndu6VvLclv1Tv3tXVJbj9XtznMZsipSOovniXT7w749EP4b4tzTeuiycXII3xQYYGp2PrjwGQi9DYglEpqCrzyJR/2KrsiRqH0RXDEPXj0XX3YiuOJBEnS1O2h5ZIc9zp98wEn9J46tuf4mf3197onWr5CEJ/Ao6/DG5cJeUkizi+yEVN6HR2Q39501oBGKfND9ccwtoPY23OgvB+utQTbX92UY/G/sK1t8KhBu+Qslf1/0FjS9v8+sz+ckKeZ4bcvjujHn8YnpvtyVuj5sttqrkonvPZPjpBzodzbTAVfpbpNsUpPN/kG6TcXW6C3GVIu4tSbltmfjB1aP58chMftrxaCOJGkj88LMZtH4Sye10mzWW3A7PtCvWR14Ahh5RxdAjqpyOYTaDiB+8AxofDBwG629uuMre8EzAD4GDmp/E3RViKTbLAHCVbyJBnJRvAmhyvQ/TrtgVuSkaGltIYu01JFaOJLHuRjT+XU7bF1cZ0vkx8AwkORTNB55fIF0mpBy1IqXnkHKzjOCRyS3Qfq6twMENbaQQsE9z7Y1dkZuioJGP0dWnNsywi0F0Flr/JHR5AvFsk7Mc4t0W6fpCw1A0F+Lu0vJzg4cn32xq7wQkOdU7MBwpv6oV7WyPlv4Oah8ieXPVBbihwyUNXTymPREnZghWVVVpdXV1zts1xSvxw5EQm9vkqIBvH1yd73ckU2uphiG+CFzdkE12qTT52eg8NPQa4EGCwxDP1tkJaXJCRKararN+VrsiN+2eagRin6d6pGECT34T8UMbPzWIdzvEu12GE5l8Y33kpgh4gBbG3EvmFvc3xilWyE27J+KC4G+ApjcUA1Ay0olIxmSUFXJTFKT8cvD9EvCDdODHIX9Sdk5W2lNVErVPkFhxEInvdyex6kw0Oj8rbRmTkT5yERkPHA4sV9VBmTinMZkkEkQ634/GFianyXu2QdwpJuFkiNb8HWofBhrGjEcmo6umQZdnEc9WWWvXFKdMXZE/BAzL0LmMyRrx9EH8+2S3iCdqoPZBNhTx5FHQUHLDCWMyLCOFXFUnA6sycS5jCl58YXJXmGYSEJmd8zim/ctZH7mIjBKRahGpXrEi9dKsxrQL7i1SbO0FIODpm/M46dDEGjS2ELVp/XktZ4VcVcepapWqVlVWZm9dXmOcJq7OEDiE5qNk/A3T7vOfJtaRWH0Wunxv9Icj0OV7kah/1elYpgU2asWYLJCKGyF4FMli7gVXD6TjbYhvZ6ejtYquPg/C75Gc3l8PugrW/jG5FK/JOzaz05gsEPEhFX9Fy69IrkEuFQWzfrzGFkJ0FhBt8kgYrXkA6XSbE7HMz8jIFbmIPA58CAwUkcUicnomzmtMoRPxIa6OBVPEAUh838LNWk29jZ1xXEauyFX1xEycxxiTBzwDQZtejQN4wT8k53HMplkfuTEFROPfofUvo5Fpm9zura3EVQ6lp9N4LXQ3SBlS8rustGnSY33kxhQAVUXXXw91E0Aa/tu6OkGnhxFP74y3J2UXgmcbtPZ+SKwG/95I2QVZ3QnetJ0VcmMKQehlqPsPEAYNJ4/F69E15yBdX8x4cyICwcOR4OEZP7fJPOtaMaYAaN2/aTzlHyABsYVo7FsnIpk8YoXcmEKgNamPixu0NrdZTN6xQm5MIQgMp/lMUQAPeLbNdRqTZ6yQt0HNmlpmvvUJX89Z6HQUUySk5BRw9+SnkSRuIIBU3IiI3eoqdvYvYDNNuOkZHrnmSbx+L/FonJ4DenDdy5fTpUcnp6OZdkxcZdD1Wah/Hg2/C+4eSMlJtpmyAUBUNeeNVlVVaXV1dc7bTdfUl2fwt+NuJVwX3nDM5XYxYLd+3Dn1RgeTGWOKgYhMV9Wqpseta2UzPP2PlxoVcYBEPME3cxax9MvvHEpljCl2Vsg3w9oV61Ied3vdrF/VwqgCY4zJMivkm2HIiCq8/uaLCWlC6bdjHwcSGWOMFfLNcsxFh9GpewW+QLKYi4C/xMe5t/0eX8DncDpjTLGyUSuboUOnMv456/94/p5X+ejlmXTt2YmjLzyMHYYOdDqaKRAaXwEo4u7mdBTTjtioFVOUNPYlJNaBd3tEAjlpT9dcDLGvkgfcfZGOf0e8NpnHtF5Lo1bsitwUFY0vRVePgtjChlUEE2iHK3GVHJO9NrUeXXkS6Bqg4cIpvgBdNRIq306OETcmDdZHboqGqqKrToPYAiCUXL9E62DdNdndizL0GhBmQxFPpklu3hB6JXvtmqJhhdwUj9inkPgOaLohQxitfSR77caXgYZSPFAHiWXZa9cUDSvkpngkVpH6n7xC/PvstevdCVL1w0tJ8jFj0mSF3BQP704t7EUZgMAB2WvXNzS5D2aj1Qv94O4Pvn2y164pGlbITdEQV0coO4fGe1H6wd0dCR6bvXZFkM4PQ9lZ4O4Nrl5QeibS5d+I2H9Bk76MjFoRkWHAbSTX1rxfVW0FKZMTGluI1j4AsXng3REp/T3i7tni811l56LeQWjtw8m9KAO/RkpGIq7SrOYUCSBl50PZ+VltxxSntAu5iLiBu4CDgcXANBF5XlU/S/fcxvwcjc5GV50CGgFiEJ2D1j8FnScg3pYnaYl/X8S/b+6CGpNlmfhcNxhYoKpfqWoEmAAcmYHzGvOzdO1fksMHiTUciYLWouuvczKWMTmXiULeE1i00feLG441IiKjRKRaRKpXrFiRgWZNMVONQ6yFD32R6bkNY4zDMlHIJcWxZvP+VXWcqlapalVlZWUGmjXFzUXqPSwByW5/tzH5JhOFfDHQe6PvewFLM3BeY1okIlByLM2LeQBKRjoRyRjHZKKQTwMGiEg/EfEBJwDPZ+C8ph1775mpnLvHZZzQaxTXnvB3Fn+x+e/90uFP4N8H8IN0SP4aOAgpOzfjeY3JZ2mPWlHVmIicD7xKcvjheFX9NO1kpt16+vaXGH/54xu2zZs88UM+mjSTu6tvoteAHq0+j4gf6XQ3Gl8CsW/A0x9xb5Gl1Mbkr4zMRlDVl1V1W1Xtr6o2ZMC0KBKO8tAVExrtfaoJJVwX5t9/m9imc4q7J+Lfy4q4KVo2rczk1LKvvk95ezwRT/Dpe/NyH8iYdsDWIzdo5GM0NAnEjQQOQ7w7ZK2tTt0riEXiKR/r1rdr1to1pj2zK/Iil1h3I7rqZKh7EGofQFeeQKLm7qy1V965A3sf/csN+57+yF/i56TLj85au8a0Z1bIi5hGP4O6x4AQyaH/ieTva+5BYwuz1u6l95/NvscOxRvwEijx06FzGRfcdTq7H7xz1to0pj2zrpUipqE3gUiqRyD8NnhOzUq7/qCfyx6+gPPvOJ31q2qo7NUFt8edlbaMKQZWyIuZeEl+KGu6Y46r4bHsKi0vobS8JOvtGNPeWddKEZPAoSSH/jel4D8413GMMW1khbyIiacPlF9BcmZkMLn1GH6ouAFx23o4xhQK61opcq6S41H/ARB+B3BD4FeIq5PTsYwxm6GgCnl9TT0L5y2lS4+OdO3Zxek47Ya4KxsWoDLGFKKCKeSPXf8Uj133NG6vm2gkxi6/GsQVEy6mpENw0z9sjDHtWEH0kb/75Ic8dsMzhOsj1K2rJxqKMuutT7j5d3c5Hc0YYxxXEFfkT4x9lnBtuNGxaDjG1JdnsH51DR06lTmUzGxMVfn0g8/5Zs4ieg7Ygp33/wUuV0FcKzSi0XnJewYShMAwxN3d6UjG/KyCKORrlq9NedztcbF+VXYLeSKRYNV3aygtDxIsy143Tt36et6Z8D6L5y9lwG792fvowXh92R/LnSn1NfVcdsjf+PqThWhCcblddOvTlVveuYaKruVOx2sVVU3u91n3HyAKeGD9/6EVY3EFhzsdz5gWFUQh3/WAHXnj35NJxBtPXPEFvHTvm71hch88P43bzrmPmtU1qCp7/WYwl9x3dsYL+uIvlnLhXmOIhKKEasMEywI8dNUE7phyPeWdO2S0rWwZP+ZxFsz8hmg4uuHYkvnLuO2c+7jqyUsdTLYZotOg7kmSSxYANCzutfYy1L834mr734VqHKLTIbEWfLsjrs5pxzXmRwXxufeUq4+jpDyIx/vT5BW3x40v6OfKI29i5lufZLzNz6u/5PqT/sGqZauJhKJEwzHef3Ya1534j4y3NfZ3d7F+VS2hhu6j+poQyxf+wPjLH894W9nyxiOTGxVxgFg0zgfPTyMeT73aYTaoJtDEGlRjm/+z9c/zUxHfiLgh/N+2Z4otQFfsj64+C117Gbp8PxI1/2zz+YxpqiAKefe+lYz7+BYOO+tgem/XE48v+UHih8UrmTZpJleOuImX7ns9o20+MfZZIvWNC1M0HGXmm5+wfNEPGWunvqaeL6q/RLXxftWxSIzJEz/MWDvZFoumLpwaVzTRbC/urEjUPYEuH4ou3wtdvgeJmjtQbbr8QBtJqj3GN001ga46HRLLQWtBa4Aw1N6NhqdkJpspegVRyAEqe3Xh/NtPZ49hu4Aq8dhPV3nhujD/vPRfREKpFoBqm6Xzv2tWXAG8fi8rFq3MWDvicqXcaAHA7c79X08ikWj0Z9taQw7fHVeTvCLCL/YeiMeb/R48rX8Z1l0HuhqIJotmzf1o7T2tPocERwCBFCePg2+ftgWLfgK6juTqko0Co3WPte2cxjRRMIX8R9WvziIWbV5oRISF85a0+byfT1vARftcyWElJ3FCr1EEOvhxe5r/8UTDUfru0KvN7TQVKPGz0747NCuCXr+XA0/eN2PtbErNmlpu+O3tHFYykuGBE7l43ytZMOtrnr1zEqP3HMMfD7qGyRM/TPnmBnDWLafSqXsFgdLkrvb+Eh9lnUq5ZNzZOcmvNbfTvFukHmrvT/ZPt4Z3Dyg5jua3jhRin7U1GC2+U2vqm/jGbK6CuNm5sc49OrFwbvOCHYvG6FjZttERX89ZyB9+dTWhhn0kVy5dzfrVtYjbhSuhJBq6Bvwlfo4aPZyyjqVtfwEp/PHB87ho7ytYv7qGaDiG1+eh9/Y9OfWa4zNy/vraEP+dOIXlC39g4B792f2QnRsNC1RVLjvkr3w9eyGxSLKL5NP353HeHn/G43Nv6GKaN3U+s96ew+i7zmzWRtctO/PgvNt467H3mD/jK/r+ojcHn7xfxv+sWpT4LvVxDYPWg2x6ZJOIQPBotK7pvYkwuvos6PYBIpt5o9u7C6Tsrw+C30bCmMwouEJ+7KUjmDtlfqPNez0+NzsMHdjmafuPXfcU4SbdMpH6CF6/lz1/swdz3ptHRddyjr10BIf8bv904qdU2asLD8+/g49ensl3Xy9n6537svP+v0gWljR9O3cxl+x7JZFwjFBtiGBpgD7b9+Tmt64mWJrsRpg7dT4L5y4hGvmp4KiCxhNE6n/qYw7Vhnn1wbc55uLD6blN893ug2VBDhvl0KqJngEQ/bj5cVdHkNa/mWj900ALN0rDkyHw682KJa5StHxMstuHCJBIjk93b42UHLVZ52qPNLYArX0AYgvAuytSepptot0GaRVyETkWuBrYHhisqtWZCPVzBg/fldOuP5Hxlz+O2+MiFomx/ZBtufLJS9p8zvkzvk55Q87r93DyVcex9U5904ncKh6vhz2P3KPZ8WgkyqJ5SynvUtamN6obTrqN9atq+LFHpL4mxNefLOSJsc/xu4Yr/iVfLGv1m4bL7eKTyXNTFnInSYc/JW8qNupeCUDZHzfvDVHraL4+O4Amr+zbwFVyPOr9RfJKP7ES8R8CwcMR8bXpfO2FRj5CV51Bcsx+HKKfofVPQZeJiKef0/EKSrpX5HOAo4GcjqU6evRhHHrGQXwzZyGdundMeyx53x16sXTBMpp2/8YiMSp7O7c412v/eoe7Ro9HVYlH48k3rP9c0uoJNqu/X8PCeYubva5IKMob/3p3QyHfalDvDd1Hm+JyuahoYxdWNolvD+j8ALr+/yA2H9w9kbLRSGDzPiFI4BA09FJDQd+IxsC3V9vzeQchFde1+efbI117JY3feKOgMXT9WKRT629SmzRvdqrqXFX9PFNhNkegxM92gwdkZELQSWOOwRdsfHXkL/Fx0Mn7OTb9f87787j93PuoW1dP/foQkVCUOe/P4y+/GZuZBja6SB2w29YM3KN/ow2RW5pa7w14qfp1fu6tKb49cHV5Alf3Gbi6vrDZRRxIjk7x7dWwNjsk/4sEoOxiW6M9gzRRA/FU+8IqRKbmPE+hy9moFREZJSLVIlK9YsWKXDXbKgOr+vPXZy+j17Y9EJcQKA1w5HnDGH3XGY5lmnjrC0TqG/fbx6NxFsz8miULlrXqHJ26d6TP9r2aDYH2BbwcfMp+jY5d99LlHH7WIZRWlOAL+hgyYncuvf8cyjqWUlIeJFAWoPtWldz85l8KaumAzSXiQjregVT8HYJHQfBEpMujuMpOdzpa+yI+Uu9OBUhhzGbOJ9LScLINTxB5A0h192GMqj7X8Jx3gD+0to+8qqpKq6uz3p3eJtFIFI/Xk5Ebjek4b/Cf+aL6y2bHS8qDXPfi/zJo7+1bdZ6F85Zw8T5XEg1HCdWFCZT46btDL25+62oCJf5N/nwsGuOL6V/hD/rYeqe+jv+5mPYjsfbPUP8ijTcAD0LZaHvjbIGITFfVqqbHN9lHrqoHZSdSfsqXq82qQ3bm608WNpv2Ho/G6bcZN1/7bNeTR7+9h/8+lRx+uN3gbdj1wB1bvSqhx+thhyHbblZ2Y1pDOlyFxldCZEryCl3DEDwCKf2909EKTsENPywWR114KJMeeIv1q2s2jO0OlPoZecUxm73zfKDEz8En77fpJ5qMSX7STSDSQveBQVwlSOf70NhiiC8GzzaIu6vTsQpSusMPjwLuACqBl0Rklqpu3kBbk1LHygr+Oetmnhj7HB9NmknHbuX8zyVHsOeI5kMU2ytVZfKTHzJh7LOs+X4tux64I6dcfRxbbNXN6Wgt0kRdcinc+ueAKOrdFSm/BvEOdDpa3hJPL/BkbrZ0MdpkH3k25HMfuckfj147kQk3PbthVUiX20VJeZBxH99CZa/83LM1sfJkiM6kUb+vlCFdJ9kGFSZtLfWRF9xaK06qr6nn0esmcsaOl3De4D/z6kNvk0hkaHU900jd+noeu+GZDUUcIBFPEKoJ8cTYZx1M1jKNftEwu7TJ4m0aSTHt35jMsT7yVoqEo4ze8wqWLlhGJJS8AXnnBQ/w8buf8qcHz3c4XfuzaN4SPF43kSaTKWPROLPfbeMCVtkW/xrE02yhQ4hAbK4TiUyRsCvyVpr85Id89833G4o4JNceefc/H7L4i6WOZPpq9reM/f2dXLj3FYwf8xirv1/jSI5s6LJlp0Zrv2wsb/vIPf1bWCDLD55BOY9jiocV8laa8eZsQjXhZsddLuHTD3I/uXXaKzMZveflvPnIZD774HMm3voCZwy6mOUL82uyVVt17dmFXX41CK+/8YdGf4mP4/50pEOpfp54tgHfYGDj8fkC4kdKTnIqlikCVshbqVvvrht2JtqYy+Wi8xYdc5pFVbn1zHsJ10U2rJESDceoWVPHQ1c9kdMs2XTFhIv55aG74/V7CZT6Ke/agT88cC6D9trO6Wgtkk53QcnIhtmJXvDtg3R5EnHn581Z0z5YH3krDT/9QCbe+uKGMd2QXL86WB5kt4N2ymmWlUtXsW7l+mbHE/EE1a+mWMq1QJV0CPKXp/5AzZpa1q+qoVvfrrjd+T0uW8SPlP8Zyv/sdBRTROyKvJU6ditnt4N23PC9iLDF1t245e2rcXtyW1yCHYItrlbYoXOONnLIobKOpfTYunveF3FjnGKFvJX+dtytTH/tp6tdVWXN92tzsh9lU6XlJfzy0N3wNunqCZT4Oeaiw3OexxjjLCvkrbBkwTJmvjWn0YgVSPZLP3P7S45k+uOD57L90G3xB32UVpTgDXgZfuaBDD/jQEfyGGOcY33krbD4i2V4fZ5my8rGojEWzPzGkUylFaXc8vY1LJ6/jOULf6Dfjn3o1K3CkSzGGGdZIW+FPtv3bLYKIYDH52HgHts4kOgnvQb0oNeA/Np2zRiTW9a10go9+nVn8PDdmu0i5At4+c0FthO6MR+mxDEAAAuNSURBVMZZVshb6fLHL+Q3FwynrGMpbq+bXQ4YxG3vX9fmxZvisTifV3/Jt58twomFywrJ59Vf8s8/PMy9lz7M59MWOB3HmLxjqx86YMqL0xl76h3EYwkS8QRde3fhb89dRq9tt3Q6Wt4Zf8VjPP2Pl5I3mhV8QR9HXTCc028Y6XQ0Y3LOVj/ME0sWLOPaE25l/epa6tbXE6oLs+SLZVz6q6uJx+JOx8sr385dzFN/f4lwXQRNKKpKuC7MM7e/zDefLnI6njF5wwp5jk26/01i0cYFW1UJ1YSY+dYch1Llp6kvTieR4s0tFo3z4fPF+4nOmKaskOfYiiWriEebFydVbVerF2aCx+fB5W7+T9TllmaLaRlTzKyQ59gev96FQGnz3evjsTiD9s7fxaCcsO//DEl5XETY99ihOU5jTP6yQp5j+x47lJ4DeuDfaChjoNTPr087gB79bCuwjXXt2YWLx52FL+AlUBYgUOrHF/Ay+p4z6dbbNuk15kc2asUBobowz9/9Ku9MeI9gWYAR5w5j32OHIiJOR8tLa39Yx5QXpwMw5PDdqeha7nAiY5zR0qgVK+QmqyLhKDNen024LswuBwyyImxMGloq5GndMRKRm4EjSO42+yXwe1W1O3YGgM8+/Jwxh91AIpFAVYlH45x+40iOHn2Y09GMaVfS7SN/HRikqjsBXwD/m34k0x5EwlEuP/R6atbUUreunvr1ISKhKOP/9zHmz/jK6XjGtCtpFXJVfU11w26zU4Be6UcqbKu+W82Tt77AfZc9wrRXZ5FIJJyO5IgZr89GU2x+EQ1HeWX8Ww4kMqb9yuRg3NOAFjeMFJFRwCiAPn36ZLDZ/DHzrU+4asRNJBIJIqEoz9/zGtsP3obrJ41xZAMKJ4VqQyjNC3kiodSurXMgkTHt1yavyEXkDRGZk+LryI2eMwaIAY+2dB5VHaeqVapaVVlZmZn0eSQei3Pt8X8nVBfesAFFqCbEZ1Pm88r4tx1Ol3u7HDAo5cSnQKmffY5JPT7cGNM2myzkqnqQqg5K8fUcgIicChwOjNQiXsZvwcyviUaar1kergvzxiPvOpDIWR0rKzjt+pPwl/hwuZLDKgOlfnbcdweGHLG7w+mMaV/SHbUyDLgM2E9Vi/rzssvtIkVPAgBub3FuGnzMRYczaO/teWX8m9Stq2efY4YwdESVbaJsTIal23F7J+AHXm+YzDJFVc9OO1UB6r/LVpSUl1BfE2p0PFDqZ/jpxbuP5sCq/gys6u90DGPatbQKuao6u89ZHnG5XFzz7J+47OC/kogniEVjuNwuho6o4oCT9nY6Xs7FY3GevXMSL977OpFQhH3+ZwgjxxxDh05lTkczpt2xmZ0ZVl8b4oNnp7F2xTp22n8Httmln9ORHPHXY2/ho0kzCNclN6z2+jxU9u7CuNm34A82XzTMGLNpWZnZaZoLlgY4cOQ+WW9n/oyvmPLSdPwBH/sfvyfd+uTPSKBvPl3E1JdnEKmPbDgWjcRY9d0a3n78fYaddoCD6Yxpf6yQFxhV5a4Lx/PK+LeIhqK4PG4e/ssTXDzuLA767X5OxwNg3kcLNoxU2VioNszH735qhdyYDLNlbAvMJ/+dy6sPvk24LkIiocQiMSKhKH8fNY71q2ucjgdAt95dkBSF3Ov3sGX/LRxIZEz7ZoW8wLw94T3CdeFmx90eF9MmzXQgUXM/rnLYdHcfj9fD8DOKdwSPMdlihbzAJItjinXLBcSVH3+dLpeLW9/9K9sPGYDX78EX8NJj6+7c8MoVdN2ys9PxjGl3rI+8wBx40j68+uA7za7KE7EEg4fv4lCq5ip7deEf/72WNSvWEglFqezVxTbOMCZL8uMSzrTaDkMHctTo4fiCPrw+D/4SH76gj8v+dQGlFaVOx2umY2UF3Xp3tSJuTBbZOPICtejzJUx9aQb+oI+9jxlCp24VTkcyxmSZjSNvZ3oP7EnvgT2djmGMyQPWtWKMMQXOCrkxxhQ4K+TGGFPgrJAbY0yBs0JujDEFzgq5McYUOCvkxhhT4KyQG2NMgbNCbowxBc4KuTHGFDgr5MYYU+DSKuQi8jcRmS0is0TkNRHZMlPBjDHGtE66V+Q3q+pOqroL8CJwVQYyGWOM2QxpFXJVXbfRt6VA7tfENQb49rNFfPDcNJZ++Z3TUYzJubSXsRWR64BTgLXAr37meaOAUQB9+vRJt1ljAKivqeeqI8cyd8oXeLxuopEYg4fvxuWPX4jX53U6njE5sckrchF5Q0TmpPg6EkBVx6hqb+BR4PyWzqOq41S1SlWrKisrM/cKTFG7a/R4Pv3gc8L1EWrX1RMJRfnolZn8+28TnY5mTM5kbIcgEekLvKSqgzb1XNshyGRCPB7niLLfEg3Hmj1W0bUDE5ePdyCVMdnT0g5B6Y5aGbDRtyOAeemcz5jNkYgniEfjKR8L1YZTHjemPUp31MqNDd0ss4FDgAszkMmYVvH6vPTftV+z4yLCrgfu6EAiY5yR7qiVY1R1UMMQxCNUdUmmghnTGhfdO4pgWQCPzw2A1++ltGMJZ996qsPJjMmdotx8OR6P8+i1T/HM7S9Tu7aObXbtx/m3n8YOQwc6Hc1spm1378/9c27l2Ttf4avZ37L9kAGMOOfXdOre0eloxuRMxm52bg6nb3beds44Xn/kXcJ1kQ3H/CV+7phyPf0G2dBIY0x+ysrNzkK0btV6Xnv4nUZFHCAaivD4Dc84lMoYY9qu6Ar5sq+W4/E171FKJJQvP/4m94GMMSZNRVfIe2zdLeW4Y5dL6L9TXwcSGWNMeoqukJd37sAhp+6Hv8TX6Lg34OPEy492KJUxxrRd0RVygAvuOoPj/ngkZR1LEZewzW79uOm1K+1GpzGmIBXlqBVjjClENmrFGGPaKSvkxhhT4KyQG2NMgbNCbowxBc4KuTHGFDgr5MYYU+AcGX4oIiuAb3PYZFfghxy2lyv2ugqLva7Cko+vq6+qNtsr05FCnmsiUp1q7GWhs9dVWOx1FZZCel3WtWKMMQXOCrkxxhS4Yink45wOkCX2ugqLva7CUjCvqyj6yI0xpj0rlityY4xpt6yQG2NMgSuaQi4iN4vIPBGZLSLPiEi72GZdRI4VkU9FJCEiBTFU6ueIyDAR+VxEFojIn53OkwkiMl5ElovIHKezZJKI9BaRt0VkbsO/wQudzpQJIhIQkY9E5OOG13WN05k2pWgKOfA6MEhVdwK+AP7X4TyZMgc4GpjsdJB0iYgbuAsYDuwAnCgiOzibKiMeAoY5HSILYsClqro9MAQ4r538fYWBA1R1Z2AXYJiIDHE4088qmkKuqq+p6o+bdU4BejmZJ1NUda6qfu50jgwZDCxQ1a9UNQJMAI50OFPaVHUysMrpHJmmqstUdUbD79cDc4GezqZKnybVNHzrbfjK61EhRVPImzgNmOR0CNNMT2DRRt8vph0UhmIgIlsBuwJTnU2SGSLiFpFZwHLgdVXN69flcTpAJonIG8AWKR4ao6rPNTxnDMmPhI/mMls6WvO62glJcSyvr4QMiEgZ8BRwkaquczpPJqhqHNil4V7aMyIySFXz9h5HuyrkqnrQzz0uIqcChwMHagENoN/U62pHFgO9N/q+F7DUoSymFUTES7KIP6qqTzudJ9NUdY2IvEPyHkfeFvKi6VoRkWHAZcAIVa1zOo9JaRowQET6iYgPOAF43uFMpgUiIsADwFxVvdXpPJkiIpU/jmoTkSBwEDDP2VQ/r2gKOXAn0AF4XURmici9TgfKBBE5SkQWA0OBl0TkVacztVXDzejzgVdJ3jj7j6p+6myq9InI48CHwEARWSwipzudKUP2Ak4GDmj4PzVLRA51OlQG9ADeFpHZJC8uXlfVFx3O9LNsir4xxhS4YroiN8aYdskKuTHGFDgr5MYYU+CskBtjTIGzQm6MMQXOCrkxxhQ4K+TGGFPg/h+SHWAEP0HFVAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "x_train, y_train = create_toy_data()\n",
    "plt.scatter(x_train[:, 0], x_train[:, 1], c=y_train)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Logistic回归"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\nabla E(\\vec w) = \\vec Z^T(\\vec y-\\vec t)\n",
    "$$\n",
    "$$\n",
    "\\mathbf H = \\mathbf Z^T \\mathbf R \\mathbf Z, \\quad \\mathbf R =(y_1(1-y_1),\\cdots,y_N(1-y_N)\n",
    "$$\n",
    "$$\n",
    "\\vec w^{new} = \\vec w^{old}-\\mathbf H^{-1} \\nabla E(\\vec w)\n",
    "$$\n",
    "$$\n",
    "\\delta \\vec w =\\mathbf H^{-1} \\nabla E(\\vec w)\n",
    "$$\n",
    "$$\n",
    "\\mathbf H \\delta \\vec w=\\nabla E(\\vec w)\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "class LogisticRegression(object):\n",
    "    \"\"\"\n",
    "    Logistic regression model\n",
    "\n",
    "    y = sigmoid(X @ w)\n",
    "    t ~ Bernoulli(t|y)\n",
    "    \"\"\"\n",
    "\n",
    "    @staticmethod 静态方法\n",
    "    def _sigmoid(a):\n",
    "        return np.tanh(a * 0.5) * 0.5 + 0.5\n",
    "\n",
    "    def fit(self, X, t, max_iter=100):\n",
    "        \"\"\"\n",
    "        maximum likelihood estimation of logistic regression model\n",
    "\n",
    "        Parameters\n",
    "        ----------\n",
    "        X : (N, D) np.ndarray\n",
    "            training data independent variable\n",
    "        t : (N,) np.ndarray\n",
    "            training data dependent variable\n",
    "            binary 0 or 1\n",
    "        max_iter : int, optional\n",
    "            maximum number of paramter update iteration (the default is 100)\n",
    "        \"\"\"\n",
    "        w = np.zeros(np.size(X, 1))\n",
    "        for _ in range(max_iter):\n",
    "            w_prev = np.copy(w)\n",
    "            y = self._sigmoid(X @ w)\n",
    "            grad = X.T @ (y - t)\n",
    "            hessian = (X.T * y * (1 - y)) @ X\n",
    "            try:\n",
    "                w -= np.linalg.solve(hessian, grad)\n",
    "            except np.linalg.LinAlgError:\n",
    "                break\n",
    "            if np.allclose(w, w_prev):\n",
    "                break\n",
    "        self.w = w\n",
    "\n",
    "    def proba(self, X):\n",
    "        \"\"\"\n",
    "        compute probability of input belonging class 1\n",
    "\n",
    "        Parameters\n",
    "        ----------\n",
    "        X : (N, D) np.ndarray\n",
    "            training data independent variable\n",
    "\n",
    "        Returns\n",
    "        -------\n",
    "        (N,) np.ndarray\n",
    "            probability of positive\n",
    "        \"\"\"\n",
    "        return self._sigmoid(X @ self.w)\n",
    "\n",
    "    def classify(self, X, threshold=0.5):\n",
    "        \"\"\"\n",
    "        classify input data\n",
    "\n",
    "        Parameters\n",
    "        ----------\n",
    "        X : (N, D) np.ndarray\n",
    "            independent variable to be classified\n",
    "        threshold : float, optional\n",
    "            threshold of binary classification (default is 0.5)\n",
    "\n",
    "        Returns\n",
    "        -------\n",
    "        (N,) np.ndarray\n",
    "            binary class for each input\n",
    "        \"\"\"\n",
    "        return (self.proba(X) > threshold).astype(np.int)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAATYAAAD8CAYAAAD9uIjPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dd5yU1dn4/881fdml7LICCijVAioWRKMmFkSxR01RfxpbYkw0scaaJ5bE5xfTfGJMNFgSE3ssCUGNYtRo7IhYAOkiCEhZlrJl6vX9456FLTOzMzuzU6/36zUvmbuc+6zAxTn3uc45oqoYY0w5cRW6AsYYk2sW2IwxZccCmzGm7FhgM8aUHQtsxpiyY4HNGFN2LLAZY3qdiNwvImtF5OMk50VE7hCRxSLyoYjs1+7cOSKyKP45J53nWWAzxuTDn4GpKc4fC4yNfy4E7gIQkTrgRuBAYBJwo4jUdvcwC2zGmF6nqq8CDSkuORn4izreAgaIyI7AMcBMVW1Q1Y3ATFIHSAA8uah0pupq63To0OGFeDTBSASvL/0fOxiLAOD3pn+Pait+tzfjuhmTK3Nmz12vqjtkU8bkow/Vhg2N6T5vLtDa7tA0VZ2WweOGAivafV8ZP5bseEoFCWxDhw7nqcefLcSjAViyvoHhIwelf33TOkYPrkv7+tbIIgBG9cvqz5UxPVYXGL882zIaNjTy0huPp/u8VlWdmMXjJMExTXE8pYrtiq5Ytjaj65d8kaoV3VHAMzbT6hhT6VYC7btxw4BVKY6nVJGBbXR9+q0vgNHVTssrk+AGsHTzuoyuN6aCTQe+FR8dPQjYpKqrgeeBo0WkNj5ocHT8WEoVGdjaZNJqawtu6WprtVlwMwZE5BHgTWA3EVkpIheIyEUiclH8kmeBpcBi4B7g+wCq2gD8FHg3/rklfiylgrxjKwaj6+tYsr6BFcvWZva+7YuGtN+3BTxjt71vM6aSqeoZ3ZxX4OIk5+4H7s/keRXdYutplzRT1mozJr8qOrC1ybRLagMJxhS3ig9smbba2mQa3KzVZkz+VHxgAye49eZAQhsLbsbkhwW2diy3zZjyYIEtznLbjCkfFtg6yUdumzGmd1lga8cGEowpDxbYOrGBBGNKnwW2JCy3zZjSZYEtgXx0ScFabcb0FgtsKdgkeWNKkwW2JNpabZbbZkzpyVlgExG3iLwvIjNyVWahWW6bMaUply22S4H5OSyvaFhumzGlJSeBTUSGAccD9+aivGJiuW3GlJ5ctdj+D7gaiCW7QEQuFJFZIjKroWFDjh6bH5bbZkxpyTqwicgJwFpVfS/Vdao6TVUnqurEurqB2T62ICy3zZjSkIsW2yHASSLyKfAocKSIPJiDcouK5bYZUzqyDmyqep2qDlPVEcDpwEuqelbWNStSlttmTPGzPLYMWG6bMaUhp4FNVV9R1RNyWWaxsQ1gjCl+1mLrIRtIMKZ4WWDrActtM6a4WWDrIcttM6Z4WWDLkg0kGFN8LLBlwSbJG1OcLLDlgE2SN6a4WGDLUr5y26zVZkz6LLDlgOW2GVNcLLDlkOW2GVMcLLDliE2SN6Z4WGDLod7ObbNJ8sakxwJbL7DcNmMKywJbjllumzGFZ4Gtl1humzGFY4GtF9gkeWMKywJbL7FJ8sYUjgW2Xma5bcbknwW2XmS5bcYUhgW2PCj23DaNrEFbZ6Ktr6Ma7HE5xhQLC2y9rJg3gFFVdOufoPFS2HofNP0BGs5Hw4t6VJ4xxcICWx4U7ST58BwIvgCEgRBoi/PZ/DNUoz2qgzHFwAJbHhXdQELrTEjY9YxA5JPMy+shjSxFN92IbjgL3fgDtPXVvD3blCcLbHlSlLltGkpxLpx+OVnQyDJovAHCH4E2Q/Rz2HoX2jI9L8835ckCWx4VXW6b/zDAn+BEDLx79OjZGWt+FOjcagxC82NoLwVX1SDa8gK6+efo1nvQyIpeeY4pHAtsBVA0Awn+g8E7DiQQP+AGfFD9A0QSBbxeEFmc+LjGIJpZ2ks6NNYCjVdB058g9A60Pg+NP0Jb3875s0xHIjJVRBaIyGIRuTbB+dtFZE78s1BEGtudi7Y7121z3pPrypvURtfXsWR9+n9hR1fvwJKmdSz5ooHRg9Pvzi7dvI5R/VK3+ETcaL8fQ/h9CM0C6Qv+IxDPkLSfkzXXEIhtTHAiBq7+uX9e6wyIrsUZMIk/hxA03Yn690fE/kr0BhFxA78HpgArgXdFZLqqzmu7RlUvb3f9D4B92xXRoqr7pPs8+10skBXL1jJ85KC0rm0LbukKeMbSGkkvZUNEwLef8ymEPt+AzT+nY3fUB4HJiCuQ7K6eC73J9qDWXgyin4FnVO6fWaKC0XAuk78nAYtVdSmAiDwKnAzMS3L9GcCNPX2YdUULwDaA2U58E6DmYnANwPl31geBY6D6/F56YlXiwxpNfq5CiQQIeMam9QHqRWRWu8+FnYobCrR/mbkyfizBc2UXYCTwUrvDgXi5b4nIV7uru7XYCqSnXdJMpdMlLTQJHIr6DwHdChJAxNt7D6s6HrYu7ZTmIuAZgnh27L3nlr/1qjoxxXlJcEyTXHs68IR2TKbcWVVXicgo4CUR+UhVlyR7mLXYCqzoctsKREQQV9/eDWoAvi+B/2jAA1LlfNw7QM11vftcsxIY3u77MGBVkmtPBx5pf0BVV8X/uxR4hY7v37qwwFZANkk+/0QEqTkPau+Cmkug3w0w4C7EM7jQVSt37wJjRWSkiPhwgleX0U0R2Q2oBd5sd6xW4sP0IlIPHELyd3OABbaCsw1gCkPcAxH/lxDvOGcAxfQqVY0AlwDPA/OBx1V1rojcIiIntbv0DOBRVW3fTd0DmCUiHwAvAz9vP5qaiL1jKxKZjJICGaV/ZDJKWoqcvwMRwGNBqoip6rPAs52O/aTT95sS3PcGsFcmz7IWWxGwDWB6TltfhI0XwIZvwsbz0ZYXCl0lUwSyDmwiMlxEXhaR+SIyV0QuzUXFKpFtAJMZbX3ZWW4pFk9Qj22Cpj+hrTMLWzFTcLlosUWAK1V1D+Ag4GIRGZeDcitKUU6SL3ZJ55k+WojamCKSdWBT1dWqOjv+6y04LwYTJt6Z1Ipuknyxi61PcnwjHd89m0qT03dsIjICJ7/EZhRnwXLb0uROkqLh3sEGESpczgKbiNQATwKXqermBOcvbJtu0dCwIVePLTuW25aBqm8Bvk4HfVB1diFqY4pITgKbOOniTwIPqepTia5R1WmqOlFVJ9bVDczFY8taKea2qYbRlmfQxqvQxqvR1ld6tUsogYOg7xXgHgZ4wL0T9L0UCRzaa880pSHrPDZx2vz3AfNV9TfZV8m0zSMtpdw21Shs+nF8fbV4MNt6h7P8+IBbMyxLIboEImvAMxLxJH9lK/5J4J+0/d5YC9r0EARfA3GDfzJUndj7U7VMUclFgu4hwNnARyIyJ37s+ngynumhkpskH5oFkSV0mdccmY+2voIEDk+rGI1thc03Q+RzEAGNor59oe+V3a6VphqBzdc79xJxDjY/DuEP0X432nu3CpKLUdH/qqqo6t6quk/8Y0EtR0pmICH0Ds6ijQm0PJ1+OVvvgshyoNXZMYsQhN6HloRvODrV4V2IfsG2oOYchPACKOOZF6Yrm3lQxPKV27aocS0fvzafVx59gyXvL+vZezHp/BK/HW1KqwjVsBOcOgQmgJCzhHd3wp+AtiY4EbXAVmFsrmiRG11fx5JeXG23cd0m7v/pA7RsacWzshWX283QsTvyw99/G18gg/dSgROSBx/v/umVoRGSLtGVaketNu5BOKOkna4VL7h69o+EKU3WYisRvbXa7l9u+htL344Rbg0TDkYINgf57JOVPHvPixk9TzxDwX9kgjMB6PO19MpwVYFnWKIz6S1d7v+KM2DQ+V7xgy/VGoim3FhgKwG9NUm+tTXEwtlLicWcd2ORoc7OVJFghLdmvJdxPaXvJVBzObiHgvQD31dgwO2IO4PBiZqL47tmdepMSH232/GJqy/0/6mT9oHXKcMzEvr9r42KVhjripaQnG8AE9ve7Vu7uIpBY1q2n4okGQjohgS+DIEv9+heAPGMQfv/BhqvAbbEjyq0PgPRT6H//3Rz/yiovRONbgBciLu2x3UxpctabCWiNzaACfTxM3zXnWjLgli7uIrIUD8ej5t9Jme0/FVuRRbTdSepEETmocn2Ie1E3AMtqFUwC2wlpKdd0lS+dfM3qOpbhS/gjGp6A14Cew3k5O8f06M65kRkPpBgdFNjELbRTdM964qWoIy7pF+sSzojYccRg/jpP67h7efmsO6z9ewybhjjjqiiekCfXFY5M67BJB7d9ICrvhA1MiXGAluJyXRGQptU062qaqo4/Otf2va9NbKosNv2BQ6Hlsc6ZX6IM6jgS7k5kTGAdUVLUrlvACOuftD/ZnDvyPbRzdHx0U37t9h0z/6UlLBSmiSfKfGMQQfcCbEGwI24B3Q4rxqD8DxnpoF3HOIqYNfZFB0LbCWqp5PkUwW3mCqfvL2YVYvXMGiXekYf6C5ol1REwN11iSuNLIPNt2yfjaARtPoCpOroPNfQFCsLbCUuV7ltLVtb+M2Ff2T95w1EwhE8Xg99a2v44X3HQL9c1jg7qhHYdBPolo4nmu5HvWMRz8juy4h8Cq3PQXSDM6PBfyTiCvRKfU1h2Du2IqDR9WjLP4g1PZZ2nhbkdpL80797ji8+XUuwOUg0HCXYHKRhTSNP/O+c4lptN/wRXSfJA4Shtfut97T1dWi8FlpfhPBsaPoLbLoSjbV0e68pHRbYCizW+hq68WJnccSWx9HGHxPbenfaK2zkaiBh1vNziISjHesWjfLx658QjWnxBDdNFoAUdGvqWzUMTXfhpJG0/f8NQXQ9tM7IYSVNoVlgKyCNNcHWO3Gy7NtWtghB8D8Q/jijsrJdty1pHFUIuMdkVJde5RkfXwWkswB4D0x9b3Q5iVcPCUPorRxUzhQLC2yFFJqTYDUKQINo6NW0i8lFl3TC4eNxuTvWxeUSdp04Grfb+WNSDK02cfeHPqcD/nZH/c5kd/9B3dzdBzSa5Fx1bipoioIFtkKSVP/7M/+tyaZLeuplxzNgUD/8fZyA4evjp3pANWf++DSg8Llt7UmfU5zJ8P6vOMsR1XwX+t/cbY6beHYCz45A5yXC/VB1fK/V1+SfjYoWkncfZ/5jF37Ef1hGRWW7AUy/2hpufPIq5vz7I1YuWs3gEYPYb8reBALbV8YtdG5be+IdB95xmd/Y93rYfBPENgICGobA8eCb1N2dpoRYYCsgcVVB3yvRLb+OH4kBLghMdf7iZijb3Dav18MBU/flgKnOtKVIJML7L33EF8vXM3TsEMYdvBuQww1gCkDcOziJv5FFEGsEz65dkn9N6bPAVmDiPwC80yD0ZjyLfj/EMzyrMnOR29a4fjO/Ou/3NG1uIdwaxhvwUju4P1fe933cVSuzql+hiQh4dy10NUwvsndsRUBc/ZDAMUjVyVkHtVzltj1865M0rt1MsDlILBYj2Bxk3coNPP3bZwh4xhbFuzZjkrHAVoayzW2LRmPMe3PhtiXDtx0PR5k988Nt3y24mWJlXdEyltUk+XbpXl6fcsDkTew5qZmtm3xoZHlGAwmq6mT5t850Xtb7vwL+Q5FEqS7G5IAFtjKV7UDC7pPG8Mk7i/H4Ivzot59RNziMPwCxWBO66Rq05hJwD942kKCRzyHyIVADvgM6zr1s+pMT1Ag638PzIPgftN//2O7spldYV7TM9bRLeuYNp9J3YA1HnrqVgUOcoAbgcgGEoOkufPiJhWaxePVt0HgFbP0zNN0NGy9AwwsA0Mjq+H6jwXZPCkLkEwi/n+2PZ0xCFtjKWDYbwNQNqeWWv1/NEae68PkTXKQh2HQVgZan43sUhJ2PtjifTT9Gm/8O4Vl0TYjFGQEOZb7FnzHpsMBW5rLZAMbr81JdmyxfLYozv9WZt7mkpfPaRlFofhiaHiVhYMMD0jejuhmTLgtsJaBlawsL3l3MyoWr0171o7OeTpKXwHF0nJfZVSDpH6MI0ELHbmgbF/iPSLtOxmTCBg+K3IsPvso/73oBj9dNLBqjbsdaLr7jfOqGpJ8tn80GMKMGH+BMOWqdDuLFWR4ojNNi2y6AiyUt/RhdtTlFiQGcTUwVqn+AeAZnXCdj0mGBrYh98vYiZtw9k3AwTDjobCD8xfJ13HX5n7nhkcsyKmt0fR1L2qV/LJq9lBf/+h8a1jSy26QxTDn7MPrXb+9Oto2SCoJUn4VWHQ+RBeDqj4bmQcvjdN3UmBTBze2sVhuYAt49EPEluMaY3LDAVsRefvR1Qq0d99aMRWOsW7GB1cvWsmMGOWptVixby+dzP+PR2/5OuLUtWK7lnWdmc/0jlzFgh/4drm9L/xBXLfjiywK5R6Kh/0L0C5yNjd2AEJA+tGpjkidHgTDim5BxnY3JlL1jK2JbNiZeEdbldtG8qSnj8kbX1xGNxnji1//cFtQAouEYzVtamHH3zI7XxwcSOk+3EgkgA36B1FwE/sOh6hRkwJ3IgN+Ca1CCgQQAH3j3zrjOxvSEBbYitvdh4/H6uzaqY9EYw3Yf1qMyN65pJOLqmvGvMeWtGbNYv6pjEEu2lLjgRfxfwVXzQ1x9zkTcgxBXf6pq/wje8XTsDHjA1R/8R/aozsZkygJbETvs61+if31/fH4v4Lx39wW8nHbFCfgD3h6VucfwIcSiMfB1fcelMeWJX01PeF+iDWASEYSqfj9jCd8Czy7gGgSB42DAr5xlmozJg5wENhGZKiILRGSxiFybizINVNUEuO6hH3L8RVMYs+9I9p28Fz+489t8+dRu1vZPoV9dDXuO2Cnp+XlvLexyLNOd5AXB5TuQZa7rkbq7kZpzEZflrFW67uKEiJwrIutEZE788+12584RkUXxzzndPSvrwQNxZjL/HpgCrATeFZHpqjov27KNE9ymnH0YU87ObEXdVM679Qwu/+YvUJ8PQh0HJ9yexBPTndy2dSWzk7wpLhnEicdU9ZJO99YBNwITcZZneC9+78Zkz8tFi20SsFhVl6pqCHgUODkH5ZpeUt2vikMO3AOXp+Nvv8fnZtKx+6a8N1mXVAkSa51JbOvviLU8jaqT8mFLG5m4bOLEMcBMVW2IB7OZwNRUN+Qi3WMosKLd95VAl76SiFwIXAiw045Dc/BYk41vXHUii3+0ng2rG5GwM0K646jBnHJp8k1Nkq22q7FN6KYfQWwLziwDL9ryBP5+PyNIJK2lxFWjzmbIsSbwjkPctdn8eCYHguFI2u9WgXoRmdXu+zRVndbue1pxAjhNRL4CLAQuV9UVSe5NGURyEdgSTQTsMu8n/kNOA9hrzwk9mxdkcqaqpopb/vA9Xn1zPhtWNTB+4khG7rUziX87O+qwbhugzQ/FN0dpm40QBg2jW+8gMOD2brukGlnubLCi8W6xhtE+X0P6fKNHP5vGtkLonXZLrQ/pUTmVzu/yZPJ+db2qTkxxPp048U/gEVUNishFwAPAkWne20EuuqIrgfbrWQ8DVuWgXJOmlq0tLJq9jC8+W5/RfSLCYQePY8Lh4xi51y6kE9QS5raF3qbzFCsAoitRdfLtknVJVRU2/xRim7avDEIEmp9GQx9l9PMAaOh9aPgObL0Xmh6AxkvRpocyLsfkXLdxQlU3qGrbxOJ7gP3TvbezXLTY3gXGishI4HPgdODMHJRr0vDsvS/y/J9exuP1EA1HGbbbTlz0m3OoGZD+BsANXzSy5IPlTDhoDLVDaomEwoSCYfr0rSJRsOvSJRVvin8/XakHEiIL4sGssyC0/gt8e6X9c2isFbb8ki6T7ltmoL79EO8eaZdlcq7bOCEiO6rq6vjXk4D58V8/D/yviLS9nzgauC7Vw7IObKoaEZFL4g93A/er6txsyzXdm/Pyx8x84D+EgxHCQWf5oOXzVnLftQ9z6d3f6fb+YGuYe6/+KwvfW0J0cA1P/bKF/vX92Lx+M6pQO6g/Z9xwGrsfMCbh/du6pP4p0PIUHeeOusG7FyJO7pqzAcyiru/aNNHKH23nmrv9GToIzyFxqzMIwVfAAlvBJIsTInILMEtVpwM/FJGTcJaFaQDOjd/bICI/xQmOALeoasqXfznJY1PVZ1V1V1Udraq35qJM070XH3yNYEvHdI1oJMqSjz5l0/ot3d7/xK+ms/C9JYSDEWKfNRLFTcOqjURCUaLhKOs/b+DuKx5g1ZI1Xe5t/+5F+pwa37zY73wkAO5BSM0PutzXpUvq2RU0QTdW/OA/pNufoaNEm0+3iWRYlsm1RHFCVX8SD2qo6nWqOl5VJ6jqEar6Sbt771fVMfHPn7p7ls08KGFbk8wldXvcNG1K3dqJRmO889zsbS29bTrNSIiEIvz7wVcTltG2bpvgxdXvRqT/rUj1BUjf65y5o66OI5sBz9guZYirCmouBHxsb20FwL0L+DPM3fNOAE0UwALg+3JmZZmSZqt7lLDxB+/OhlVvEo10bPG4XS4G71Lf4ZhGN6EtD0LwbWIxN58vHxdfV2171827vplwfZ+O98VirPk0dS7athVAPKPAM2r7veGP0ab7IfoZuPpC1amoZ/cu6R8SOBL1jHI2fIltAt+B4D8IkcymjYmrGq35Pmz9A07rLQr4wf8lJ+iZimGBrYQdc97hvDdzDs1bWomEIoiA1+/lG9ec3GEGgcZa43lmTkqGS2DITm9y4U1V/P76BJPp281IcHtdjJqwS9I6jK7egVlLFvLKqwsIVPuZcNg4qqoDaGQRuvlnQLyrHNsEzQ/jD5xEyDexa3DzjICa7t8LdkcCh6Ge3SD0mpPu4ZsInt1tN6wKY4GthPUb2JcfP3YFLz/yX+a/vYi6IQOYfNZhjNyz427yGnwNYltpn5Lh9cHocS0MH9PKisXOFlTiku2tNp8PiYTxBXxMPjNZN0554vYZvPTWeyDg3xLmsdv+zvduP5cxYx5hW1DbXhFomY6/z2kEI8tz9z+iE/EMAc/Xe618U/wqPrA1bWqmeUsLdTvW4naX3ivHmgHVnPi9Yzjxe8ckvyj6Cc6CkB0pMGxMkFXLqtlp7BB2O2AMgWofb06fRaMXdt9nV776w2MZMKh/l3sB5r+1iNeffgdpCRHe2QfNzgjntKv+wm1PrEicFScCUWdAK50ZCcVEY81Ot1rqEE/mi3ya/KnYwNaytZW/3PQ4c99YgNvtwuP3cvo1X2X/KWW4GKJrGM7L+Y4tKFVhw2of3oCHr11xImP3GwnAcd8+atseCTsMqyeZN/45i1B8VNb7WYjwzlV4GlqIxWI0Nw2kuk+COcoaA3ctAYaUzCR5VYWWv0HzEzh/ZSKodw/oew3i6tPd7aYASq+JkiP3XPsgc9/4hEgoQrAlRFNjE3+9+XGWfth7XaRCkcBkEHeHHNpoBDZt8LDoQ6cbOqJT9zWdbfuioa4jkJE6J29t7brJOMG0PT8Ejkfiu145uW0lMEk++F9ofhwnZaTV+W/4Y9jymwJXzCRTkYGtYc1Glry/jEio42hiOBhm5l8TpzaUinUrNzD3zYU0rNneWhJXP6T/z8A9gljUCWqfvN+H3127C16fj3Nu/iZeX9fG++j6upTb9k06dl98VduDl/czp/UWiynDxh2J9LsO3PHBCamBPl9Hqs/qUk7RB7emP9M1R04h/D4a6z5f0ORfRXZFG9dtweP1dMnhUoUNqzLfpq4YhFrD3HPtgyx8d7Hzs4UjTDh8T869+Ru4PW7EMxJ37W+Ihrcw9/XFfPzeUg46uYaDTzyA+mGpW2cr2u1u1d6EI8azx7/GMv+tRYRaQri9btTrYcolx8ZX+J2ADLgDRZEk81BLYt22pMt+qTMoY4toFp2KDGw7jhxEJNw1293tdbPr/qMS3FH8nrh9BgvfXdxhetVHr8zlX/e9xPHfnbLtOre3LxMO35cJh6ded61Nqj1JRVxc+IuzWThrKR/9dz7V/aqYdNx+bOwb6bACSLKg1l6pDSRs47a9UYtRRXZFq2oCHH3u4fjbdaNcbheBPn6OyuFKtfkSiylvz5jVpQUaCob5z9/ezMkzkndJhV0njua0y05g6vmTqRtSm/FS4olmJBQVz8jEx91DEKnIv0JFr2J/V4779mTO+snXGb77UGoH9+egE/bjuod+yIAdEm0dV9w0FuvyvrBNsCXFJPM0tQ0kpHrflkgGixQW90BC9UV0HQjxQs3lhaiNSUNFdkXBWYts/yl7l0V6h9vjZthuO7Hik887HBeBsfvmpmudqkua8Pokq+12pxi7pOIdiw74FbQ+DZFPwTMCAqcgnp5tgWh6X8W22MrNGdeegr+Pf9tUKrfXjb+Pn9OuPDGnz8mk1dY2ST5dxdwlFc8wpOYHyIBfO/+1oFbUKrbFVm5G7DmcGx6+jH8/8hqfL1rNiPE7c8Tph1A7OPGsgZ7ItNXWpvNS4t0pxlabKS0W2MpI/bA6vvmj3t0gbHR9HUuSpH8kvD7DLmlb+ocFN5MN64qaHuntgQRjsmEttiLWtKmZD16ZSzgUYfzBu1E/NP3uXC5EIzFeefx1/vvk24SDYfY9ai+mnndkjwcSrEtq8sUCW5H68NX53H/dQ4hLiMWUp/5vBkefdwTHf/uovNXh/hseZu7rnxBqdfYy+M9jb/DhK/O44VEnzSHZjIREetolNaYnrCtahFq2tnD/9Q8RCoYJtoQIB8OEgxFm/vkVls9bmZc6rF76BR//d3tQA4iEo2zesIVZz89Ja5J8ImWT22aKmgW2IjT39YVIgrXhwqEI7zw7Oy91+HTuClyurlOhgi0hFs1eCnQ/Sb6zTGcktLHgZjJlga0IRSPRxPt0qhKJJJ5hkGu1gwdAgsDm8Xm6vOur1Nw2U7wssBWhcQfvSjTaNYD5Aj72n5KfTUl2nTiamgHVuDq1HN0eN4d8ddK27/nokoK12kxmLLAVob61NXzzRyfj9XudJYdE8Ff5mDh1n22r3PY2l0u4YtpFjNxrZzw+N16/h/qhA7nkdxckXCq8N7ukba02C24mXTYqWqQO+eokxuw3ilkvzCHcGmbvw8Yzcs/hed1tqXZwf66893ts2biVcDBC7eD+CZ/flsz86NUAAA64SURBVP6RySgpZDYjwUZJTSYssBWxwTvX5zW9I5m+tTXdXlPJk+RN8bGuqMkpG0gwxcBabGVAVXnvhQ959Yk3CQdDTDxmX7582kH4ApntpJ6tfEySd3LbFlmrzaRkga0MPPz/P82s594nGN8Kb/WStcx6fg5X3f/9DjvC50NvT5JvY11Sk4oFthKz9KPP+OcfnmfVktXsMLyeQ0+ZxDvPvNdhWfBQMMyaT9cx56WP2f/o/KSHdGYDCaaQ7B1bCVk0exl3fG8aC95dzJaGJpZ+sJyHbn0KTZDMG2wOMu+thfmvJJnntrWlf1hum8kVC2wl5Kk7nukwdxMgGo4SDXfduNjtdSfMN8unfOS2GZOIBbYSsmrR6oTHVZ39Ddpzu10cfNIBeahVYhW/AYwpKAtsJaTfwMQb8/qrfNQPHYi/ykeg2k91/z5857azGbhTbZ5r2FFPu6SZsuBmOstq8EBEfgmcCISAJcB5qtqYi4qZrqaeP5knfjWdYGto2zFfwMeUcw7j2Asms3rpWsLBMMN23Qm3p3j+zcp43bYv1tlAgslKtn/6ZwJ7qurewELguuyrZJI5+OSJHPfdKQSq/fgCXvxVPo44/WCmnj8ZEWGn0YPZZdywogpqNkneFEJWLTZVfaHd17eAr2VXHZOKiDDl7K9w5BmHsHnDVmpqq/H6ij9jxzaAMfmWy3/azweey2F5JS0WU9Z8uo4NqzfmvGy3x03t4P4lEdTasw1gTL50+zdDRF4EhiQ4dYOq/iN+zQ1ABHgoRTkXAhcC7LTj0B5VtlQseHcxf/6fR2ltChKLxRg8YhAX3nY29cPyuxlLMbENYEw+ddtiU9WjVHXPBJ+2oHYOcALw/6kmShXdVs40VZ2oqhPr6gbm7icoMg1rNnLX5X9m0/ot8f0KIny+aDW3f/duotFYoavXq75Yvo5n7pnJ3+98jmUfr0h4jeW2mXzIdlR0KnANcJiqNuemSqVLVfnrzU90SaLVmNK8tZUF7y5h3EH5/csXjcZoWL2RPjVVVA/o02vPefWJN3ny9meIRWLEYjFeefQNDjxhP06/5qvb1nCzSfImX7J9SXMn4Admxv/wvqWqF2VdqxL15j9nbdvopIuYsnn95rzWZ/aLH/HIz58mEgwTjUbZbdJYzr3ldKr7VeX0OZs3bOHJ22d0nK/aGuKdGbOZNHUfRu+zfdVfmyRv8iGrwQNVHaOqw1V1n/inYoMawL/ue4lYku5mNBpl5J47560uy+eu4C83PUZTYxPBlhCRUJRP3l7EtB/9JefPmvvGAlyurn+UQsEws2Z+lPAeW7fN9KbiSXgqA5s3bE16bp8j9mLwiPy1IF588NUOLShw5pUu/3gF61ZuyOmz3B531zldOIc83q5/xCy3zfQ2C2w5NHy3nRIer+5XxTm3fDPj8lSdlJHVy9aSYlwmofWrGhLe4/a62bQut13iPQ/dHU3QUvV4PUyaul/S+2wDGNNbLLDl0CmXHd9l1VpfwMuZP/4a7gQbIKeyfP5KfnLybdx29h384lu/439O+jnL5yYeaUxkt4ljcHu7LjIZDkfYcXSi7J2e69O3ivNuPQOf35kN4fV78fo9TL3gSIbvnjjY52uSvCkeIjJVRBaIyGIRuTbB+StEZJ6IfCgi/xaRXdqdi4rInPhnenfPKq0MzyI3aq+dufyei5hx1wusXLia+uF1nHDhFHY7YEzaZUSjMT5+bT733/AI4eD20dVgS4jffv8ebn3mOqpqun/5f+SZh/LGP96lZUvrtj1K/QEfR5z15ZwPHgBMOHw8tz5zPR/8Zy7hUIQ9D92dgTumnoRvG8BUDhFxA78HpgArgXdFZLqqzmt32fvARFVtFpHvAb8A2ro6Laq6T7rPs8CWY7vsMYyL7zi/R/duWL2R279z97bt7jqLRZXZMz/ikFMmJbi7o34D+3L9w5fy7L3/Zt6bC6juX81RZ32Ficf03oq61QP6cPDJmS+VZJPkK8IkYLGqLgUQkUeBk4FtgU1VX253/VvAWT19mAW2InLvtQ/SuHYTsVji92mhYIhNG7akXd6AQf058/pTc1W9XmG5bcUrHIpk8qqgXkRmtfs+TVWntfs+FGj/LmUlcGCK8i6g4xTNQLz8CPBzVf17qspYYCsSjes28/miNUmDGoC/ys/ofUbkr1J5Yrltxcnv8WQygr1eVSemOJ9op++Ef9hF5CxgInBYu8M7q+oqERkFvCQiH6nqkmQPs8GDIhEOhnG5ku/y7gt4GTFuOLvuPyqPtcovG0goayuB4e2+DwNWdb5IRI4CbgBOUtVg23FVXRX/71LgFWDfVA+zwFYk6ofWUd2/65QnEehbV83JlxzLxb87b9v0pHJjG8CUvXeBsSIyUkR8wOlAh9FNEdkX+CNOUFvb7nitiPjjv64HDqHdu7lErCtaAE2bmvnbr//J7Bc/RGPK+EN245tXf5Vzf3oGf7j0fqLRKJFQFH+VjwGD+nP1AxenNRJaDjIeSLB120qCqkZE5BLgecAN3K+qc0XkFmCWqk4HfgnUAH+L/wP+maqeBOwB/FFEYjiNsZ93Gk3tQjJN/MyFvfacoE89/mzen1sMYjHl1tNvZ+1n64lGnDQMl9tF39oabvr71TRvauaNf7zD+lUb2W3iaPafsjdef353dC+ktoGEdINbW2DLZGmj1kj5DyTUBca/1807r25l8vd01/HDsn5eLlmLLc/mv72Ihi8atwU1gFg0RmtzK7NnfsCXTpzI8RdOKWANC8ty20wu2Du2PPti2VoinZY1Agg2h1i1ZE0BalScbJK8yYYFtjzbsHpjwgUnvQEPO+V4qlOpsknyJlsW2PIoGonx9oz3Ep5Thf2m9N6sgFIzur7OJsmbHrPAlkeN6zYRDUcTnquqCeAPVM4gQbost830hAW2PKru14eoJl6Ism5IYXdtL0aW22Z6ygJbHgWq/RxwzD54/R0Ho30BH1PPO7JAtSp+tgGMyZQFtjw7/dpT2H/KPnj8HvxVPqpqApxy2XFMOHxcoatWlPIxkOBMkrdWWzmxPLY88/o8fOumr/P1q05ka2MTdUMGOEtrm6RskrzJlLXYCqSqJsAOwwZaUMuA5baZdFlgMyXBcttMJiywmZJiuW0mHRbYTMmwDWBMuiywmZJiuW0mHRbYTEmy3DaTigU2U3Ist810xwKbKUm9PUm+jQW30mSBzZQ0y20ziVhgMyXLcttMMhbYTMmz3DbTmQU2U9Ist80kYoHNlLye5rZlylptpcMCmykbNpBg2lhgM2XBcttMezkJbCJylYhofPt5YwrCcttMm6wDm4gMB6YAn2VfHWOyZwMJJhctttuBqwHNQVnGZMUmyRvIMrCJyEnA56r6QRrXXigis0RkVkPDhmwea0y3bJJ8Zes2sInIiyLycYLPycANwE/SeZCqTlPViao6sa5uYLb1NiapfOW2WauteHUb2FT1KFXds/MHWAqMBD4QkU+BYcBsERnSu1U2pnuW21bZetwVVdWPVHWQqo5Q1RHASmA/VV2Ts9oZkyXLbatMlsdmypZNkq9cOQts8Zbb+lyVZ0wu9HZum02SL07WYjMVwXLbKosFNlP2LLet8lhgMxXDctsqhwU2UxFsknxlscBmKoZNkq8cFthMxbHctvJngc1UFMttqwwW2ExFsty28maBzVQc2wCm/FlgMxXJJsmXNwtspqLZQEJ5ssBmKpbltpUvC2ymolluW3mywGYMNpCQDyIyVUQWiMhiEbk2wXm/iDwWP/+2iIxod+66+PEFInJMd8+ywGYqnk2S730i4gZ+DxwLjAPOEJFxnS67ANioqmNwNom6LX7vOOB0YDwwFfhDvLykLLAZE2eT5HvVJGCxqi5V1RDwKHByp2tOBh6I//oJYLKISPz4o6oaVNVlwOJ4eUl5clr1NH0898P1u44ftrwXiq4HSmWxy1KqK5RWfUuprtA79d0l2wI+nvvh87uOH5buJugBEZnV7vs0VZ3W7vtQYEW77yuBAzuVse0aVY2IyCZgYPz4W53uHZqqMgUJbKraszew3RCRWao6sTfKzrVSqiuUVn1Lqa5QvPVV1ak5LE4SPSLNa9K5twPrihpj8mElMLzd92HAqmTXiIgH6A80pHlvBxbYjDH58C4wVkRGiogPZzBgeqdrpgPnxH/9NeAlVdX48dPjo6YjgbHAO6keVpCuaC+a1v0lRaOU6gqlVd9SqiuUXn0zFn9ndgnwPOAG7lfVuSJyCzBLVacD9wF/FZHFOC210+P3zhWRx4F5QAS4WFWjqZ4nTkA0xpjyYV1RY0zZscBmjCk7ZRnYROQqEVERSTcHpyBE5Jci8omIfCgiT4vIgELXqbPupsEUExEZLiIvi8h8EZkrIpcWuk7dERG3iLwvIjMKXZdyUnaBTUSGA1OAzwpdlzTMBPZU1b2BhcB1Ba5PB2lOgykmEeBKVd0DOAi4uMjrC3ApML/QlSg3ZRfYcOaYXU03CXzFQFVfUNVI/OtbOPk5xSSdaTBFQ1VXq+rs+K+34ASMlBnqhSQiw4DjgXsLXZdyU1aBTUROAj5X1Q8KXZceOB94rtCV6CTRNJiiDRTtxVeG2Bd4u7A1Sen/cP4RjhW6IuWm5PLYRORFYEiCUzcA1wNH57dGqaWqr6r+I37NDTjdqIfyWbc0ZDyVpRiISA3wJHCZqm4udH0SEZETgLWq+p6IHF7o+pSbkgtsqnpUouMishcwEvjAWRCAYcBsEZmkqmvyWMUOktW3jYicA5wATNbiSyrMeCpLoYmIFyeoPaSqTxW6PikcApwkIscBAaCfiDyoqmcVuF5loWwTdEXkU2CiqhbtKg8iMhX4DXCYqhbdYl3x+XoLgcnA5zjTYs5U1bkFrVgS8SVuHgAaVPWyQtcnXfEW21WqekKh61IuyuodWwm6E+gLzBSROSJyd6Er1F58YKNtGsx84PFiDWpxhwBnA0fG/3/OibeITIUp2xabMaZyWYvNGFN2LLAZY8qOBTZjTNmxwGaMKTsW2IwxZccCmzGm7FhgM8aUnf8H7LqinqoNZioAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "x_train, y_train = create_toy_data()\n",
    "x1_test, x2_test = np.meshgrid(np.linspace(-5, 5, 100), np.linspace(-5, 5, 100))\n",
    "x_test = np.array([x1_test, x2_test]).reshape(2, -1).T\n",
    "\n",
    "feature = PolynomialFeature(degree=1)\n",
    "X_train = feature.transform(x_train)\n",
    "X_test = feature.transform(x_test)\n",
    "\n",
    "model = LogisticRegression()\n",
    "model.fit(X_train, y_train)\n",
    "y = model.proba(X_test)\n",
    "\n",
    "plt.scatter(x_train[:, 0], x_train[:, 1], c=y_train)\n",
    "plt.contourf(x1_test, x2_test, y.reshape(100, 100), np.linspace(0, 1, 5), alpha=0.2)\n",
    "plt.colorbar()\n",
    "plt.xlim(-5, 5)\n",
    "plt.ylim(-5, 5)\n",
    "plt.gca().set_aspect('equal', adjustable='box')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Softmax回归"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "class LabelTransformer(object):\n",
    "    \"\"\"label encoder decoder\n",
    "    \n",
    "    Attr\n",
    "    ====\n",
    "    n_classes: int\n",
    "        number of classes\n",
    "    \"\"\"\n",
    "    def __init__(self, n_classes=None):\n",
    "        self.n_classes = n_classes\n",
    "\n",
    "    @property\n",
    "    def n_classes(self):\n",
    "        return self.__n_classes\n",
    "\n",
    "    @n_classes.setter\n",
    "    def n_classes(self, K):\n",
    "        self.__n_classes = K\n",
    "        self.__encoder = None if K is None else np.eye(K)\n",
    "\n",
    "    @property\n",
    "    def encoder(self):\n",
    "        return self.__encoder\n",
    "\n",
    "    def encode(self, class_indices):\n",
    "        \"\"\"encode class index into one-of-k code\n",
    "        \n",
    "        Params\n",
    "        ======\n",
    "        class_indices: ndarray with shape (N,)\n",
    "            non-negative class index\n",
    "            elements must be integer in [0, n_classes)\n",
    "        \n",
    "        Returns\n",
    "        =======\n",
    "        output: ndarrray with shape (N, K)\n",
    "            one-of-k encoding of input\n",
    "        \"\"\"\n",
    "        if self.n_classes is None:\n",
    "            self.n_classes = np.max(class_indices) + 1\n",
    "\n",
    "        return self.encoder[class_indices]\n",
    "\n",
    "    def decode(self, onehot):\n",
    "        \"\"\" decode one-of-k code into class index\n",
    "        \n",
    "        Params\n",
    "        ======\n",
    "        onehot: ndarray with shape (N, K)\n",
    "            one-of-k code\n",
    "        \n",
    "        Returns\n",
    "        =======\n",
    "        output: ndarray with shape (N,)\n",
    "            class index\n",
    "        \"\"\"\n",
    "        return np.argmax(onehot, axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "class SoftmaxRegression(object):\n",
    "    \"\"\"\n",
    "    Softmax regression model\n",
    "    aka\n",
    "    multinomial logistic regression,\n",
    "    multiclass logistic regression,\n",
    "    maximum entropy classifier.\n",
    "\n",
    "    y = softmax(X @ W)\n",
    "    t ~ Categorical(t|y)\n",
    "    \"\"\"\n",
    "\n",
    "    @staticmethod\n",
    "    def _softmax(a):\n",
    "        a_max = np.max(a, axis=-1, keepdims=True)\n",
    "        exp_a = np.exp(a - a_max)\n",
    "        return exp_a / np.sum(exp_a, axis=-1, keepdims=True)\n",
    "\n",
    "    def fit(self, X, t, max_iter=100, learning_rate=0.1):\n",
    "        \"\"\"\n",
    "        maximum likelihood estimation of the parameter\n",
    "\n",
    "        Parameters\n",
    "        ----------\n",
    "        X : (N, D) np.ndarray\n",
    "            training independent variable\n",
    "        t : (N,) or (N, K) np.ndarray\n",
    "            training dependent variable\n",
    "            in class index or one-of-k encoding\n",
    "        max_iter : int, optional\n",
    "            maximum number of iteration (the default is 100)\n",
    "        learning_rate : float, optional\n",
    "            learning rate of gradient descent (the default is 0.1)\n",
    "        \"\"\"\n",
    "        if t.ndim == 1:\n",
    "            t = LabelTransformer().encode(t)\n",
    "        self.n_classes = np.size(t, 1)\n",
    "        W = np.zeros((np.size(X, 1), self.n_classes))\n",
    "        for _ in range(max_iter):\n",
    "            W_prev = np.copy(W)\n",
    "            y = self._softmax(X @ W)\n",
    "            grad = X.T @ (y - t)\n",
    "            W -= learning_rate * grad\n",
    "            if np.allclose(W, W_prev):\n",
    "                break\n",
    "        self.W = W\n",
    "\n",
    "    def proba(self, X:np.ndarray):\n",
    "        \"\"\"\n",
    "        compute probability of input belonging each class\n",
    "\n",
    "        Parameters\n",
    "        ----------\n",
    "        X : (N, D) np.ndarray\n",
    "            independent variable\n",
    "\n",
    "        Returns\n",
    "        -------\n",
    "        (N, K) np.ndarray\n",
    "            probability of each class\n",
    "        \"\"\"\n",
    "        return self._softmax(X @ self.W)\n",
    "\n",
    "    def classify(self, X:np.ndarray):\n",
    "        \"\"\"\n",
    "        classify input data\n",
    "\n",
    "        Parameters\n",
    "        ----------\n",
    "        X : (N, D) np.ndarray\n",
    "            independent variable to be classified\n",
    "\n",
    "        Returns\n",
    "        -------\n",
    "        (N,) np.ndarray\n",
    "            class index for each input\n",
    "        \"\"\"\n",
    "        return np.argmax(self.proba(X), axis=-1)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQgAAAD8CAYAAACLgjpEAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dd3xb5b3H8c9PR8sjjmNn7x1ISEhCCGGPMMIqtKUtpaVs6KC3dLHay20vl7aUttABpUAppVCg7D3C3iGD7L1jx5lOvDWO9Nw/juzYlmTLkWzL9u/9evkVWzo655HtfP2cZ4oxBqWUSsTV2QVQSmUvDQilVFIaEEqppDQglFJJaUAopZLSgFBKJZVyQIjIgyKyS0SWN3qsSETmisi62L992qeYSqnO0JYaxEPAnGaP3Qi8ZYwZB7wV+1op1U1IWwZKichI4CVjzGGxr9cAJxljykRkEPCuMWZCexRUKdXx3Gm+foAxpgwgFhL9kx0oIlcDVwPk5uUcMX7C6DQvrZRKZvGiFXuMMf3SPU+6AZEyY8x9wH0A0444zLz98X866tJK9ThF/klbMnGedHsxdsZuLYj9uyv9IimlskW6AfECcEns80uA59M8n1Iqi7Slm/Mx4BNggoiUiMgVwG+A00RkHXBa7GulVDeRchuEMebrSZ6anaGyKKWyjI6kVEolpQGhlEpKA0IplZQGhFIqKQ0IpVRSGhBKqaQ0IJRSSWlAKKWS0oBQSiWlAaGUSkoDQimVlAaEUiopDQilVFIaEEqppDQglFJJaUAopZLSgFBKJaUBoZRKSgNCKZWUBoRSKqmMBISI/FBEVojIchF5TET8mTivUqpzpR0QIjIE+C9gRmzPTgu4sKXXBGw73csqpTpApm4x3ECOiLiBXGB7SwfbYZv15eUZurRSqr2kHRDGmFLgd8BWoAyoMMa80fw4EblaRBaIyIKKfbXsLC3XkFAqy2XiFqMPcB4wChgM5InIN5sfZ4y5zxgzwxgzY0BxP1z2KA0JpbJcJm4xTgU2GWN2G2PCwDPAMa29aFxBkYaEUlkuEwGxFZglIrkiIjhb8a1K5YUaEkplt0y0QcwDngIWActi57wv1ddrSCiVvVLevLclxpj/Af7nYF8/rqCIdZWws3QTAGOLijJRLKVUmrJmJKXWJJTKPlkTEKAhoVS2yaqAAA0JpbJJ1gUEHAgJpVTnysqAUEplh6wOCL3NUKpzZW1AaFuEUp0vawMCNCSU6mxZHRCgIaFUZ8r6gAANCaU6S5cICNCQUKozdJmAAA0JpTpalwoI0JBQqiN1uYAADQmlOkqXDAjQkFCqI3TZgAANCaXaW5cOCNCQUKo9dfmAAA0JpdpLtwgI0JBQqj10m4AADQmlMi1Tm/cWishTIrJaRFaJyNGZOO/B0MVmlMqcjKxqDfwReM0Yc4GIeHH251RKdXFpB4SIFAAnAJcCGGNCQCjd86ZrZ6lzi6FL6Ct18DJxizEa2A38Q0Q+F5EHRCSv+UGNN+8tL9+bgcsmp20RSmVGJgLCDUwH/mqMmQbUADc2P6jx5r1FRcUZuGzLNCSUSl8mAqIEKIltwQfONnzTM3DetGlIKJWeTOzNuQPYJiITYg/NBlame95M0ZBQ6uBlqhfj+8CjsR6MjcBlGTpvRujen0odnExt3rsYmJGJc7UXDQml2q5bjaRsjd5uKNU2PSogQENCqbbocQEBGhJKpapHBgRoSCiVih4bEKAhoVRrenRAgIaEUi3p8QEBGhJKJaMBEaMhoVQ8DYhGNCSUakoDohkNCaUO0IBIQENCKYcGRBK6tqVSGhBKqRZoQKRAbzNUT6UB0YJxBUXs3tIb0JBQPZMGRCt6ckiY8GrM/hswe76KKb8KU/cqxpjOLpbqQBoQKeiJIWHsjVDxC7DXATZE90LNw1D7RGcXTXUgDYgU9biQqHmc+O1NglD3PM7WJ6on0IBogx4VEpHNiR8XgUhm9jUxJoiJ7MWYaEbOpzIvU4vW9hjjCopYtwX6jahgfXl5913b0hoM0T3xj5sIuArTOrUxYaj5OwTeAQTEi8n7FuI/Na3zqszLWA1CRKzYzlovZeqc2apH1CRyvwZ4mz3oA/8ZiCsnvXPXPACBd4EwEAJTDdUPYIIL0juvyrhM3mL8AFiVwfNlte4eEuI5FAquB2sQzl/5HMj5AuRdmtZ5TTQQC4fm7RghqHsyrXOrzMvILYaIDAXOBm4DfpSJc3YF3f12Q7zTwTvduSXAjYik/FoTqYC6f0HwMxA3+E5xaiWmiqR/lxLd0qhOlakaxF3A9UDS1qaO3Ly3I3W1moSJVmFqn8dU3en8G61q9TUinraFgwlCxU8h8J5z+xDdD3UvQsWtIIUgVuIXusemfA3VMdIOCBE5B9hljFnY0nEdvXlvR+oqIWHsMtj3Pah9DIIfOP/uuxZj74g/NrIbU/MwpuKXmJp/YyL7Ur9Q8EOIVgGRRg+GnTEVkc2Q+03A1/Q14oPciw7iXan2lIkaxLHAF0RkM/A4cIqIPJKB83YpXSIkau4DU8OB+/9YA2HN/U0OM/Ym2H8d1L0E4SVQ9xzs/z7GLkntOvZaIJjgCQORTeCZDO5BBx52DYCCWxH3iIN4U6o9ZWLz3puMMUONMSOBC4G3jTHfTLtkXVDWh0R4GdB8qLRxQqCx6nvB1AF27AEbTC3UPJjadVyDie8BIXZrkQcVN4K9+cDj0XKoezy1c6sOpQOlMiy7QyJZm7Sn4TNjorHh1QmEl6d2Gf8pTsNkE5YzfiKyFUzz2kUYQssxdmlq51cdJqMBYYx51xhzTibP2RU1Doms4j+B+JBwxx6vJzQOjCbEn9JlxNULev8K3GMAy/nwTIGC/wN7AwdqJo1fZEFkW0rnVx1HR1L2JLmXgb0V7C3OkGljwD0Cci9tOEREMP6TY6Mcw41e7AX/aSlfStzDofAOTLQOxIWI0yhp3KMgvJS4kDARsIYc9FtT7UMDogcRVw6m96/BXg+RErCGgntsfBdm3mUQ2QXhlc5fdmODdyrkXnhQ12zCPwcCrzjnbOABzwTEPaztb0q1Kw2IdrSz1GmDyKYBVCICnnHOR9JjfND7FqdNILLdCZLwQth3DSa63/lLn3cp4j0i5esaE2vodPWB3rdB9X1grwY84D/Jqd2orKMB0U7GFRSxrhJ2lm4CMhsSxt7kzKh0j0Ks9htTIu4h4B6CqX0Gap+koesyUgqVd2AKfoZ4J7dcVmOg9nEIvOjUGsQPuV9HCn/VsPhMWwZhqY6lAdGOmodEvYMNCxOthMpbwS4BcYEJY/yzIe/qdvtPZowNdU8TP64hBLWPgvc3LZ+g7kmoe+HA60011DyMkTykSeOoykbazdnO6pfP372ld/rdn1V3xcYPBA+MUwi8C4E3M1TaBKLVzdoLGolsb/Glxhioe574cAlCna5M1RVoQHSAcQVFjCtwag1+y7n3b2tImGhVbBxCpNkzQQi04wx7V36CMQ0x1hCMMcnXqTSBBGMeYqLZNkZEJaK3GJ3Ab40jEFnXJCRave0wAZLmuanNXOHiTw7WSLCbz+T3gBTA3q8DIYz7UMi/uulwafGDqzdEE8zjsIa3Y5lVpmgNooNt27gLcEKi/gNSqFG4+jp/zeNY4J2RsfIZexOm5j+Yuucwkd1Q+SunW7QJcdaJCC+mYV6HvQoqbsZEDrwPEYmNsWg+7NoLeZdkrMyq/WgNogM5jZblbNu4i2Gj+zc8nkqNQkQw+ddC5e04A5iigAdceZDz1bTLZoxx5loE5sbOb0HNv3HmbjS/rbGcnozmj5swBF6FvG8cKLf/eIzkQN1jztgKazjkfdNZkEZlPQ2IDlYfEs3V1yTqrS9fFx8S3qmYwjsg8DJEysBzWGwJuF7pF8xeHWvsrJ/pmaRhsuG5RJVPG+xNcY+Kbwb4MlfLUR1HAyJL+a1xrC9vOmlqbFER4h4K+ddk/oLBD0k8RTsRN/GzQmOPu8dkrkyq02kbRBZrcztFWiT2kSL3ROImdYnXGUrdihZ7PlRW0RpEF3FQPR8xxhjmbtzAUyuXUxEIMr64mMumTmd049f7Tmx2i9ESAf/JYI+E4JtgQuCZBHlXIlaf5OWwtzlrTdirAA/GfxLkXpr+Ktmq3WgNogs52BrFE8uX88Ci+eyorqbODrNk5w5ufOt1tu7f33CMeMZBzrkknerdRBgCbyD5lyHFjyJ9n0R6/8K5/UnCRPZBxU2NukvDzozRyl+n9B5U59CA6CT13Z0HK9WQCNo2z6xeTsBu2uMQikR5bPnSpgdbw0jctpCAqUm1qI7A605Nowkb7LUYe0vbzqU6jAZEJ6gfVdkRIbGrpgZJ0LYQNYb1+w6sLm7sMqi+h5Z7L+q5wXt02wob2ZT43FLfZaqykQZEJ6kPiXT5rXHsq6vlnvnz+NuC+cwr2UY0eqAWUJSTSyRJg+Cg/AJCts2i7dsp2/MSJm68QyIesPrGbkdaZ4xxbi+sESS8fTGRWM1FZSNtpOziFm4v5TcfrWRQ/m5sY3hz03oG5xfwu9Pm4HFb5Hk9HDFoCPO3l2BHD2xb4nNbTB04kIuffQoR4eIx65kzLJKkH8PtjOS0+oF3FvhORlytLz9nAp9A7f0QrcG5dWkeVB7wTNSFYrKYBkQnaz6qsi3sSITff/IhoUiELRUHaiTG7OexFUu5eMpU7l+4gEU7Sok0CofePh+XTp3OvQs/Ixhrm/hk1wBOHrSRHHfzWoQHCv+EuAe0qWwmvAqq/0jTXhG3s3GOqQbxgG825PXIBdC7jLQDQkSGAQ8DA3HG/95njPljuuftCZINvU7Vhn3lRBPcPmzaV8iC7aX0zc1j7sb1BCNN/9MHIzaVwSCNX7q0vB+f7x3AtOKdjULCB/6zmoSDsTdCzeMQ2ejsAJ77NcQzKb5wtU8T32VqO+FQ9ADiKmjz+1UdLxM1CBv4sTFmkYj0AhaKyFxjzMoMnLvbSyck3C5X0j6HysBg3tkYZeO+QoYUNN3q0IWLjfvLm9xygPDbpUdxVP8dfHNCJcN69wXfbMR7WMMRJrwOKv6bhv/40XKoWIvJ/xHin9m0ANH43bqcy3gguhc0ILqETGycU2aMWRT7vApnh29dnrgNDrbBcnSfIvI98RvU+C03c8aMa9gotbSy6bJ0UaKMKizCazX98RuExeVDCef+EOn1gybhAEDNP0i4K3ft3+NHRnoOIeGvl4mAa1D84yorZbQXQ0RGAtOAeQme65ab93YEE63G1DyO2f9jTMX/YkKLAWeG589POIl8jxe/243XsvC6LI4fPpJjh4/gpBGj8FlOJbG0srjhY1tFMROK+3LssBH43QcqkX63m+OGj2RMshGa9sbEj0d3Q8UNmNCyA4/5L3D222zCB7lfTqmBU2UHydSYeBHJB94DbjPGPNPSsZMPO9w8859XMnLd7qJ+hmfz2wwTrYb9P3J2yG7Yp8IHeRciOedjjGFndTWr9+4mYNsc1q8/A/LyWVC2ncpggKU7dzJ/ewmhSAS3ywUIP551DEcPH0FdaD4l+5azfu9uxvTazPEDN+B25TtdmP5z4ta5NPuugcjuFt6FFwquR7zTnePtUmeZvMgG52nJhdxvITmnp/8NUy0q8k9aaIxJewptRnoxRMQDPA082lo4qMSStkUEXoVoBU03sQlCzWPM230o9yxaRk0ojMFw7LDhjCzsww3PP4NtohhjiBrDkYOHMLKwDzluD5P692f1nt1sKPkto/wLGGtZjO0fICBRNlf6GVuwB2oedQYv5X+7aSFzLoDqB0k+6zMENQ9BLCCwtzTdLcvUQs0/MKAh0UWkfYshzp+ZvwOrjDF/SL9IPVfCtojQAhJNoIpg8cqa19gfCBCORrCjUT7etpWfvz2X6nCIgG0TjEQIR6MsLCtjeO9Cevm83PDm62zc+SyDPZ8hhIEAAH7j/Cqsryx0rhd4GxOpaHpR36mQ+yUSbszbULBGoyLrHk1Q9qCzeIzqEjLRBnEscDFwiogsjn2clYHzKgBX4vaASDTMnkDTkYnhaJRwk54JRzBi89La1dw9fx7haJQ5Q9fit+JHTTYOiTWVxazfu7bJ8yKC5H4Fiv8FkmSRGldho0ImuR2JVjjL6ausl/YthjHmQ9q0kED7CgfDBOtC5PXO7R4bsuScC6HPaVqtd7GzroCSmtS7CvfW1eKKfT9yrXDS4+pDIhiNsKCsFnGXJ1j+zoPJ/YpzK9KkXD7nNqSeNSDxPAtXHyTZStkqq3SbuRihQJhHbn2KH5/8C24+6zZ+ds6vWfp+85WYux7xTIT8y50eAckFvOAeyVt7Lkz5h+ezLMYX9W34+rNdgwlFk786GHGxtLw/dy8IE4lGWF9e3vDRwH825H4VJAdnhesc5+vGC8bkXkzCBWtzv4HqGrpNQDz8i/8w/7XPsUM2dijC/p0VPHjzo2xa3vW3lBf/6VD0EBTcAoV3IoW/5+wJR+Nzu5tU3XyW08XptSzc4vxo/W43Y/oU8Y0pUxtGXT61eQL7Q34CEecY2zhDYG0jhKIuPtgxjDuWziRqDNXBoQnXoHBuN74IRQ9Dn/ug6GEk94tNam3imwm9fuisgI3LmcuR/13Ef0p7f8tUhmSsm7MtMt3NWVlezX+f+2vCwab3tSIw5cRJXPO7b2XsWh1hXWV5SqMqt1dV8sjSxSzftYvefh9fPnQSJ44Yxe6aGt7ctJ79gQCT+w1kR00V727eRGUwSGXIuSXIscLMHryZqcU72RnI49WtY9gbzCEcdWEbCwCPy8W/vvgVcjwH2joCkfh1MlX2yapuzs62b+d+LLcVFxDGwK6tezqpVAdvXEER61IYej24VwHXH3sCu2uqeWvTRlbu3oXXsjhqyDAumjyVkG3zw9dfobSqMm5Idl3Ew0vbxvHStgOraQtRxvcuxyWGzdX9OGboyCbhAE1X365fAk9DovvqFgExYHg/InZ8673LcjFqSvfewenzHdv51QfvEY5EMMAbG9bTLzePO08/k49KtrKjpjqlNaIOK9rHTyd/hEec81guF+6Cn7b4msbrZGpIdE/dIiD8eT5mf+N43nnsQ4J1Tr+7iODxe5hzafe9341Go/z+Y2e6dz0D7Kqt4RvPPUWh399sQlZiBZ4o/zPtIzzSbABU9R0Yzz2IVZz4hcSHxEdbt/D0qpVUBAIcPmAgX588hX55eQf7FlUn6xYBAXDud06neHARcx9+l+r9NYydOorzvn8mfYe2/pdtb9k+PnlhAft3V3DoUeOZevIkLLfVAaVuWWszPDft30ednXw8wf5AoMXzC5Dn8XLdETaehF3CBoIfQO75LZ6nPiTunf8Zb27a0BBY72zeyKelJfz5zLMpzs1t8RwqO3WbgBARjj3/SI49/8g2vW7lp+u47yf/JBKJEglHWPjGUuY+/C4/uv87eP2prPDcPlKZBu5xWU0WgmmLHLebB879InleLwReJlpjJ+jSCseGebfOREfyzyULCEUOTC+PGEPADvPsqpVceYTurNUVdZtuzoMRiUT5x88fIxQIEwk7f/WCtUF2bNzF+0990smla30a+LDevZvMxkzELfFL1rpdLq45YiZ5Xi8iwsbqIYQi8S0VNl7wHp5SWUsqKrBitZDG08vtaJSlu5KsDaGyXo8OiO3rdmCH46vooWCY+a8t7oQStY2I8ONZx7Z4jEtc/OG0M/nujJmcOmo0wwp6gzHcPX8eV7zwDJ+Xbee+JTv5cMcQ6uwDt1V1tsWK8r7gnpJSWfrk5DRp76ifWg6Q7/HFD7RSXUKPDgi3zw0J/nICnXp70VxLy+MfOXQYNx53QsKahNeyOHzgIEYXF3PG2PGUVFVSVlWJbQzhaIS9dXX8+sP32LhvH39ZNYM/rZzBgj0DWLy3P39bPY3/W3x0i20cjfXNy+Ow/gNiU8oP2F3Tn5NGntxB2weqTOs2bRAHY+DIfvTu14vdJXubrM/oy/Fy3JdmdV7BGkmlLeLoocM56ktDWb5rJy+vW8OyXTvxWhZnjBnHBROdVaG2Vexn47592M0GxoUjEfxuNyGET3YN5ZNdB3bH8rvdrd7CNHb9Mcdz17yPWVi2Hbe4cFsurpl+JBP7OWtaprN9oOocPTogRIRr/nApd13zN0LBMCZiMCbK9NOmMPPMqZ1dvAZDXbks3LyV4r755BYk7g1wuVxMGTiIKQMTL+e2p7YWt0jc5OsoUOTPIVJX22RxW59lMWfMWPbW1dI3Ny+liW+5Xi83H38SVcEA1aEQ/fPysZrVKHSgVdfSLYZapytiR1j5yVoq91YzdupIBozs19lFAiAaNTx710u8//SnhIbmEa0LMuuc6Xzt+vNxtbEbdl9dHVe++GzcuAiPy8WXDp1Ib5+fR5ctIRSJIgIFXh/7A3WIuOjj93PdrGOY1L9tS9+nQodutw8dap1Bltti8vGHdnYx4rz16Pt88MxnhIM2oUAIQjbzXllEbu88zvvenNZP0EifnBzOGD2WNzdtaKgpuIAcj4ezxx1Cb7+fOWPHUxkM8LO35lJWXeUsemsi7Kqt4ZfvvcNfzjyH/vn5GX2PWqPIbj26kTLbvfXI+04wNBIO2Lz7xIdNVpGORqNs2r+PsqqqFs931RFHcuW0GQztVUCh388po8Zw5+ln0dvvLCJruVzsqK5mb6CO5qMrItEIr61fG3/SDGrckKm9HtlBaxBZrKairukDPg8Ew4Rqw6xfuIFxM8aycHspd376MeFohKgxDMrvxc3Hn8jA/PgVn0SE08eO4/Sx4+Keq7entjbh47YxbK9uOYAyQWsU2UVrEFls+KEHthfxlcbCwud0v97744fZuqec33z0PlWhIAHbJhSJsLViPz97ay7RaJStq0p4598fMv+1xXE1kWTGFhUlHJ3psywmt0MbREu0RtH5NCCy2NnXnNZkMb/mIfHYh58Rbfaf2QDV4RC33/Yof7jqrzz3l1d47FdPc/NZt1Gydnur1xzUq4Bjhg3HZx1oBHWL0Mvr45RRY9J+T21Vv1iNjqPoHBoQWWrj0i3cf/0jcd2L9SERNVH22cG4cQ3g9Mps2FxGOOCsrhWsDVFXGeC+nz6MMQZjDJ/v2M7zq1exYHtJXMhcd9QxfOvwaQzp1Yu+ObnMGTueP5xxVtzaEAAB22Zj+V72Berinss0DYmOl6l9MeYAfwQs4AFjzG8ycd6eyhjD3296lGBtsv0nIBoxHDtuDJvWLicQaTra0bajWJsr415TVV7NpnWl/HHDYnbUVGNHorgtpxvzN6fOoTDWWOlyuThn/CGcM/6QFsv59MrlPLFiGS4R7GiU6YMG86Ojj2vT4Kq2OhASB7pHtY2i/WRiXwwLuBs4E5gIfF1EJqZ73p5s55Y91FYm/4vs9ro568rZzJl0CP3y8vA0Gozksyz6lQSwqhOsXC3Ck1vWUFJZScC2sU2UgG2zo7qav3zWtslpH23dzBMrlhGMRKiz7dj+G9v587yOmeTW+NZDaxTtJxNRPxNYb4zZCCAijwPnAV1yd287bLPknRWsW7SRosF9mHXODAqKMtv33xpxSfxmuDF5vXP5zi3fwh6ex9bKCr4xeQrry8v5bHspPsvC77ZYYUeI/uAw3LsD9HmzFF+ZEzb+XB+LKnZjm/h2i/nbS9lZXc2AFMc5PL1yRZORl+DM3Py0dBs1oZAzjbyDOCGxTmsS7SATATEEaLx0dAlwVPODRORq4GqAwYOyc/PvQE2Q311xN3u37ydYG8Tjc/PqA2/z/buvZPTkjlu6rv+wYgr7945bT9Pr93Lud04nf3Qxv3j3bXaGaoj2cROORvnapMms3rObJTt3EBVABHtADru/MoqhT2zGW2lzxa+/wU3rP0163Qc/X8BNx5+UUhn3JVmMxhLp8IAADYn2kolGyiRLETV7wJj7jDEzjDEzioqSL2HWmd7457vs2rq34d4/HLQJ1gb5x88fS/oXvT2ICFf99mLyeufiy/Ph9lh4/V4mzBzLsecfxf+9/w7bqyoJ2ja1ser9EyuWsXhHWdzOWuKx6PPNw/nf529g7LRRzBg0OOl1F+9Ifd2Gw/oPSPiD91oWfTtp9aj62w3tFs2cTNQgSoBhjb4eCrTen5aFFryxBDsUP725ak8V5WX7ERdU7Klm0Kj++POab22fWUPGDuS2l29m6XsrqdhTydipIxkxaRglFRWxYdDNZmVGownT3ghE+ubQK3abdMW0GXxcknivEI+V+t+LiyYfzoLtJQRsu2HUpdeyuHr6kbhcndc5pgOtMisTATEfGCcio4BS4ELgogyct8MlW4cyEo3y0C2Ps3VVCZbbImJHOfOKU5hzefsuiOv1e5hxRtMVnarDISxxAc79v2efTbiP82NMtPicW4TxxQd21eqbl8dxQ4fzccnWJsd7XC5OGz025bIN6tWLP805hydXLWfFrl30z8vngomT2mVC18HS6eXpy8TenLaIXAu8jtPN+aAxZkXaJesEx31xJi/+9XVCgQM9AOISPF43W1Zsww5HGvbeeP0f7zBgRD+mzZ7coWUcXVjUUHtwRSBqOSFBsZchBQVsr6pqssq1x7I4t1l35XdnHsXO2hq2VlQgAlFjmNSvPxcdltrqUfX65efz3SOzY92MZLRGkR6d7t1IxI5w30//xcpP12KiUcTlIiffT111oGHNysZGHjac7/3pcuY+/C6fv7kMb46Xk756DEefdyQuV/ttHPzWxg3cu/AzQrG9MCyvRcHgXvzu1Dm8vH4Nr6xbS23Y5rD+/bli2hEM610Ydw5jDOvL91JWVcWIwj6MKIw/pjvqKdPLMzXdWwOikWAgzJ1X/ZUdm3YRCoSxPBaWZWFMNG7XLoCigYVYbot9u/Zjh5wA8fm9TD99Chff8pV2LevavXt5ce1q9tXVMqJvERccPz3hSEeVXH1YdMeQ0PUg2sE7//6Aso07G8IgEo4QCUeQBLUBy7LoM7CQkjXbG8IBIBgIseD1xZx5+ey4PTlCAWcx3OUfrqKwXwHHfXkWQ8YOPKiyji8u5sdHOwvWrqss13A4CLozWOs0IBqZ98rnCWsKLsuFeAQ79pzlsfDn+snJ9zfs5NVYNGKY9+oizrpydsNcimBdiDsu+wt7S/cRrAvhslx88lPcbiYAABDTSURBVOICLr7lKxxxempLy6vM06HbLdPJWo243Ym/HS6Xi0t/eSGHnzSJIeMHcdJXj+HnT/yQgaMGJOz5iNgR3njoHR77zbMN4yfef+pTdm8rbwiUaCRKKBDm0dueJpyga7WtWlr5WrVOZ4wmpjUInKr/s398mR2b4v+TiUDfIUVMP3Uy009t2mNxwpdn8f6TnxCx4xsww0Gb+a98zlFnTmPM1FEsenMJ4WCC+RHApy8uYOz00Qwa1fJu3uCsLfn+ls3UhsNMGzSICcV9U1r5WqVGu0ab0oAAHrjxEdbMX08k0nQkgS/Hi8fv4arfXpzwdX2HFvHtP1zC/Tf8i7qq+KHHoUCYhW8uY8zUUeTk5yQ8R6AmyDN3vYwxhn5Di/nOXZdRNDBxj8L80lJu//h9MBCORnhm9UqOGjKUHx19bENIqPQ17hoFevQQ7h5/i7G7ZC9r5q+Pa3twuYSx00dz28s3M7CFVa4PmTmWr/7kPHw58XMPxCW4Pc4tyIlfPSbhMeC0T4QCYco27eLu/3ow4bDuoG3zu08+IBSJEIo63ZvBiM1n20uYV1rShnes2qonD+Hu8QGxa8serATrF0SjhlAghMfbeiVr8gmHYqLx/6ndHouZZ04HYMoJh3LS147B7XXjz/NhJRjWHI1EKS/bR8nasrjnVuzeTfwum86CLe9s2tjwtbZFtI+eurJVjw+IgaP6E0mwvZzlsRgxcViCV8TL7ZXDZbddhNfnwZfrw+v34PG5OeuqUxk63tnIRkQ479ozufWFG7n4lq8wcHTiIckuy0XN/pq4x1vat6b+ufrNfjUk2ldPWiuzx7dBFA/uw+TjJ7Lsg1UNjYgi4PG6OfnCljfGbezwkyZy2ys3s/T9ldRWBZgwY0xDODTWu28vps2ezJ7t+3jp3jfiGi7tsM2ISUPjXjepb+LGR5/lZvbIA2tFaltEx+gpQ7h7fA0C4NJbv8YpFx1HXu9cLI/FhJlj+ek/vkefAb3bdiIDi99eznN/foXbL/kzv/zyHWxYsiXhocd/6SgK+xXg8R3IaK/fy9nXnJ6wQdPrtrjxuBPwWW58loVbXHgtixOHj+TIIdm5vkZP0Z1vO3SodYYYY/j1N/9E2YadTbo9fTk+fv6fH1I8qE/ca+qq63jvyU9Z8u4KevXJ4+QLj+PQWcn3rACoCgb5eNtWasJhpg0cxKg+8eddV1mu3Z2dIJuGbutcjCyzZWUJd13zt7iRlZbHYvZFx3P+98/s0PJoSHSObJkMlqmA0FuMDNmzfV/CORuRcISdW3Z3eHnGFRRpY2Un6G69Hd2+kfKj5+bzyv1zqdhTRb+hxXzxB2cz5YRDKV2/g6d+/yLrF2/Cn+vj+AtmcfZVpyZdNKY1wyYMSjii0uvzMGbqyDTfxcEZV1DEOh1d2Wm6w2Swbl2DeO8/H/Pk715g384KopEoO7fs5sGbHuXTlxby+yvucUZPhiPUVNTy9iMf8NAtTxz0tfoP68uUEybh9R+YVWlZLvz5fo4578hMvB3VBXX1LtFuGxDRqOGlv82N25MyFAzzzB9fjhs5GQqGWfreCsp37D/oa15664WcffVpFA/uQ35hHqMOH0nfIUXcc91DfPD0POxw+pOyVNfTlW87uu0tRqguRKA68c5UNRW1SUY+utmxeXfSuRCtsdwuTvvWiZz2rRP55/88weK3lzc0Wpas2c7CuUv4r3uuarLa1OrP1vPiX19n55bdDBjRj3O/cwaHzEx9bUjVtXS1yWDdtgbhzfHiz0+88nRurxwsK76twQ7b9B+W/pL8pevK+PzNZU16NEKBEFtWlrDy4zUNjy3/cDX3/vAhNi3bSm1lHZuWbeXeHz7E8g9Xp10Glb0a1yiyfWewbhsQLpdw5pWn4PM3nSDl9Xs4//tn4vY1rTx5fB4OnTWevkPST/N1Czdh4rcGIVgbZM2C9Q1fP33nS4SajaQMBcM8fddLaZehMe3NyG7ZPBksrYAQkTtEZLWILBWRZ0Ukq1Y+PfnC4/jidWdRUNwLgL5Dirnkl1/j2PNnct29VzNi0jBEBK/Pw7HnH8kVv8rMav15hbm4EvSGuL1uevU5sLXdrm174o4B4nbUSofOz+gasrWdIt02iLnATbGl728HbgJuSL9YmSEinHDB0ZxwwdEYYxqWfwMYMXEoN/zzWiKRKC6XNHnuYETsCLtLysktyGHKSZN44vbn4o5xuVzMPGt6w9e9+uRRubc67rjcXrnYYRu3JzNNRDo/o2vJpnaKtH4DjTFvNPryU+CC9IrTfpIFQKJp12312Wuf85/bnycSiRKxI0w4YgxX33ExD93yBHXVAUScTXkuv+0iCvsVNLxuzuWzee7Pr8b1tARqA1x/6v9y9jWnMfui49Mun+p6smUyWCZ7MS4Hkg4k6Aqb96YiUBNkxcdrsMMRJh49jl1b9/Lv/3u6yWY7Kz5Zw+aV2/jeny7HZbmI2FGGHzIEq9malyd+9WhCgRCvPfg2wUAIE3HaLepX037xnjfoXVwQt7tWS6r31/DJiwvZuWUXow4bwYw5U/HFxmboknRdU2fWKFqdiyEibwKJ1mb/mTHm+dgxPwNmAF8yKUzu6KpzMVZ8tIYHbngELAED0UiEASP6UbqujETv2uNz84N7r2l1Z/BgXYjrT/1lwhW1B48ZwM+f+FFK5StdX8Yfrrw3tgNYGF+Ol7zeuVz/8PcpKMpvuM3QkOjaUpkU1mFzMYwxpxpjDkvwUR8OlwDnAN9IJRy6qtqqOu6/4RGCgRDBmiDB2iDhoE3puh0JwwGchWuf/ePLrZ47WBdKeo6KPVUpl/Ffv3iSuupAwxoTwboQFbureOHuV4EDDZaqa+vIhsx0ezHm4DRKfsEYU5uZImWnZe+tQpJ8txJN0qqXaPm45vIL88jJ9yd8bmSKq1rVVQcoXb8j7vFIJMKSd1emdA7VdXTUEO50W+j+AvQC5orIYhG5NwNlykqhYLihjaAxg2lYmDaRguL8pM/Vc7mEL//onCbzOEQEX46XL1w7J6XyJepWrde8fNrl2T10RNdoWgFhjBlrjBlmjJka+/h2pgqWbSYdOz7h4Cef38tlt17IsEPiG169fi9zLj8lpfPPnDONq393CWOmjqSwf2+mnDiJn/zjuwybMDil1/v8Hg45alzcCFGPz82sLxy4FdVxEd1Te4VEt52LkWlFA/tw5hWzee3Bt7FDNtGowZfjZersyUyYOY4FbyyhdF0Z0UgUBLxeD2dfcxqzzjki5WtMnDWOia2sKNWSi2+5gDuvvpf9e6owkSggjJo8nLOumN3kOB0X0T0l6u1IlwZEG8y5/BQOnTWeea8swg7ZTD91ChOOHMPd//UgaxducMIBwAAC42eMTnsAVlsUFPfiv5/8CWsXbGBvaTlDJwxmxMT4BXBV91UfEpmiAdFGIyYObfKfbm/ZPtYt2thkh29wejDe/Nf7GRu+nSqXS3Q2aA/XfGewdHTbyVodZf/O/QmHRBtj2J3BORVKdQYNiDQNHD0w4UIwlsdi7PRRnVAipTJHAyJNeQU5nHThcU2mlbtcThflqd88oRNL1jrtyVCt0TaIDDj/2jkMGN6Xtx79gOr9NRwycyznfucMCvu3ceOdDlTfk6HzM1RLNCAyQEQ45rwju9zitNrdqVqjtxhKqaQ0IJS2RaikNCB6OB16rVqiAaE0JFRSGhAK0LUiVGIaEEqppDQgVBN6m6Ea04BQDbQtQjWnAaGa0LYI1ZgGhFIqKQ0IpVRSGhBKqaQyEhAi8hMRMSLSNxPnU0plh7QDQkSGAacBW9MvjsoG4wqKtCdDAZmpQdwJXA8J1oRXXZaGhIL0d9b6AlBqjFmSwrFXi8gCEVlQXr43ncsqpTpIqwvGtLR5L3AzcHoqFzLG3AfcB87mvW0oo1Kqk7QaEMaYUxM9LiKTgVHAktjeD0OBRSIy0xgTv0mk6pJ0Sbqe7aBvMYwxy4wx/Y0xI40xI4ESYLqGQ/ehQ6+VjoNQLdKh1z1bxhatjdUilFLdiNYglFJJaUAopZLSgFAp0YbKnkkDQrVKezN6Lg0IlRLtzeiZNCCUUklpQKg20duMnkUDQqVM2yJ6Hg0I1SbaFtGziDEdP7FSRKqANR1+4eT6Ans6uxCNZFt5IPvKpOVp2QRjTK90T5KxodZttMYYM6OTrh1HRBZoeVqWbWXS8rRMRBZk4jx6i6GUSkoDQimVVGcFxH2ddN1ktDyty7YyaXlalpHydEojpVKqa9BbDKVUUhoQSqmkOiQgROQXIlIqIotjH2clOW6OiKwRkfUicmM7lucOEVktIktF5FkRKUxy3GYRWRYrc0a6jZqdv8X3KyI+EXki9vw8ERmZ6TI0utYwEXlHRFaJyAoR+UGCY04SkYpGP8db2qs8ja7Z4s9AHH+KfY+Wisj0dizLhEbvfbGIVIrIdc2OadfvkYg8KCK7RGR5o8eKRGSuiKyL/dsnyWsviR2zTkQuSemCxph2/wB+AfyklWMsYAMwGvACS4CJ7VSe0wF37PPbgduTHLcZ6NtOZWj1/QLfBe6NfX4h8EQ7/owG4Sw6DNALWJugPCcBL3XE70yqPwPgLOBVQIBZwLwOKpcF7ABGdOT3CDgBmA4sb/TYb4EbY5/fmOj3GSgCNsb+7RP7vE9r18umW4yZwHpjzEZjTAh4HDivPS5kjHnDGGPHvvwUZ8n+jpbK+z0P+Gfs86eA2RLbYyDTjDFlxphFsc+rgFXAkPa4VoadBzxsHJ8ChSIyqAOuOxvYYIzZ0gHXamCMeR8ob/Zw49+TfwLnJ3jpGcBcY0y5MWYfMBeY09r1OjIgro1VAR9MUgUaAmxr9HUJHfMLejnOX6BEDPCGiCwUkaszfN1U3m/DMbFAqwCKM1yOOLFbmWnAvARPHy0iS0TkVRGZ1N5lofWfQWf93lwIPJbkuY7+Hg0wxpSBE/RAoo1MDur7lLGh1q3swPVX4FacH/atwO9x/mM2OUWC1x50H2xL5THGPB875meADTya5DTHGmO2i0h/YK6IrI4leCak8n4z+j1JhYjkA08D1xljKps9vQinSl0da0d6DhjXnuWh9Z9BZ3yPvMAXgJsSPN0Z36NUHNT3KZPL3ifcgas5EbkfeCnBUyXAsEZfDwW2t1d5Yo005wCzTewmLcE5tsf+3SUiz+LcFmQqIFJ5v/XHlIiIG+hNfPUyY0TEgxMOjxpjnmn+fOPAMMa8IiL3iEhfY0y7TVJK4WeQ0d+bFJ0JLDLG7Gz+RGd8j4CdIjLIGFMWu71KNB+/BKd9pN5Q4N3WTtxRvRiN7wm/CCxPcNh8YJyIjIol9IXAC+1UnjnADcAXjDG1SY7JE5Fe9Z/jNGwmKvfBSuX9vgDUtzZfALydLMzSFWvb+DuwyhjzhyTHDKxvAxGRmTi/P+22E3OKP4MXgG/FejNmARX11e129HWS3F509PcopvHvySXA8wmOeR04XUT6xG7xT4891rIOavH9F7AMWBp7M4Nijw8GXml03Fk4recbcG4F2qs863HuxxbHPu5tXh6c3oUlsY8V7VGeRO8X+F+c4ALwA0/GyvsZMLodvyfH4VQ5lzb6vpwFfBv4duyYa2PfiyU4jbvHtPPvTcKfQbMyCXB37Hu4DJjRzmXKxfkP37vRYx32PcIJpjIgjFMruAKnXeotYF3s36LYsTOABxq99vLY79J64LJUrqdDrZVSSWVTN6dSKstoQCilktKAUEolpQGhlEpKA0IplZQGhFIqKQ0IpVRS/w/uZ5tlopp4iQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "x_train, y_train = create_toy_data(add_class=True)\n",
    "x1, x2 = np.meshgrid(np.linspace(-5, 10, 100), np.linspace(-5, 10, 100))\n",
    "x = np.array([x1, x2]).reshape(2, -1).T\n",
    "\n",
    "feature = PolynomialFeature(1)\n",
    "X_train = feature.transform(x_train)\n",
    "X = feature.transform(x)\n",
    "\n",
    "model = SoftmaxRegression()\n",
    "model.fit(X_train, y_train, max_iter=1000, learning_rate=0.01)\n",
    "y = model.classify(X)\n",
    "\n",
    "plt.scatter(x_train[:, 0], x_train[:, 1], c=y_train)\n",
    "plt.contourf(x1, x2, y.reshape(100, 100), alpha=0.2, levels=np.array([0., 0.5, 1.5, 2.]))\n",
    "plt.xlim(-5, 10)\n",
    "plt.ylim(-5, 10)\n",
    "plt.gca().set_aspect('equal', adjustable='box')\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
